text:
# Self-Organisation to Criticality in a System without Conservation Law

## I Introduction

The idea of self-organised criticality (SOC) was introduced by Bak, Tang and Wiesenfeld (BTW)  as a possible paradigm for the widespread occurrence in nature of scale free phenomena. It refers to the intrinsic tendency of extended, non-equilibrium systems to spontaneously self-organise into a dynamical critical state. In general, SOC systems are driven externally at a very slow rate and relax with bursts of activity (avalanches) on a very fast (almost instantaneous) time scale. The standard signature of SOC is a power law distribution of avalanche sizes and in this sense is the system said to be critical. Typical physical realizations of this phenomena includes, among others, earthquakes, forest fires and biological evolution (for a review, see e.g., Ref. ).

A number of simple lattice models have been developed to test the applicability of SOC to a variety of complex interacting dynamical systems . In general these models reach a stationary critical state after a sufficiently long transient time. Not much attention though has been paid to the self-organisation process and studies have mainly concentrated on the properties of different models at stationarity. This is partly justified by the fact that some of the most studied models, such as the BTW  and the Zhang  models, display a very simple transient time behaviour. Indeed in these models the relaxation time does not scale with system size (time is defined as the input energy per site in the system) . On the other hand, there are also models with a more complex behaviour (see, e.g., Ref. ).

In this paper we investigate the approach to the stationary state in the so-called Olami-Feder-Christensen (OFC) model for earthquake dynamics . This model has in recent years attracted a considerable deal of attentions especially because it has been proposed as an example of a system displaying SOC behaviour even with a non-conservative dynamics . One of the most important question in this field is indeed whether a conserving local dynamics is a necessary condition for SOC . For example, it is well known that the BTW model is subcritical if dissipation is introduced . The presence of criticality in the non-conservative OFC model is still debated . Recent numerical investigations, though, have shown that the model on a square lattice displays scaling behaviour, up to lattice sizes presently accessible by computer simulations .

The present investigation complements previous analysis of the OFC model which were based on the study of the probability distribution for earthquake sizes. It provides further numerical support in favour of criticality in the non-conservative regime. Indeed, we will show that the model displays a non trivial transient time behaviour: the relaxation time scales with system size and it is controlled by a dynamical critical exponent $`z`$ that appears to be universal over a range of dissipation levels of the local dynamics. Moreover we will establish the presence of long range spatial correlations in the system. In so doing, we will be able to gain some insight into the mechanisms behind criticality in non-conserving systems, mechanisms that are very different from those at work in systems with a conservation law.

The plan of the paper is as follows. In section II we describe the model and briefly summarise previous foundings relevant to our investigation. In section III we define the quantities of interest and present the results of our numerical study. Finally, in section IV we discuss our main conclusions.

## II The model

The OFC model is a coupled map lattice model, where to each site $`(i,j)`$ of a square lattice of linear size $`L`$ is associated a real variable $`F_{ij}`$. In the initial state, at time $`t=0`$, the values of the $`F_{ij}`$ are chosen randomly in the uniform interval $`(0,F_c)`$. Subsequently the variables evolve according to the following two-steps dynamical rules: (i) if all sites in the system are stable (i.e., $`F_{ij}<F_c`$), they increase simultaneously and uniformly at a constant rate

$$\frac{F_{ij}(t)}{t}=v;$$
(1)
(ii) as soon as one of them reaches the threshold value $`F_c`$, the uniform driving is stopped and an “earthquake” starts:

$$F_{ij}F_c\{\begin{array}{c}F_{ij}0\hfill \\ F_{nn}F_{nn}+\alpha F_{ij}\hfill \end{array}$$
(2)
where “$`nn`$” denotes the set of nearest neighbour sites of $`(i,j)`$ and $`\alpha `$ is a parameter that controls the level of conservation of the dynamics ($`\alpha =1/4`$ corresponds to the conservative case). The “toppling” rule (2) can possibly produce a chain reaction, which ends when there are no more unstable sites in the system. At that point, the uniform growth (1) starts again. In the following we will assume, without loss of generality, a unit growth rate, i.e. $`v=1`$. A crucial point in the description of the model is the choice of boundary conditions and, in accordance with previous investigations, we will consider open boundary conditions. These conditions imply that sites close to the boundaries topple according to (2) but have a smaller coordination number.

There is a clear separation of time scales in the system: earthquakes occur instantaneously on the slow time scale of the driving. The time in the system is therefore set by the slow time variable $`t`$. By construction, moreover, the time coincides with the input energy per site in the system. We will make use of this latter observation when we will try to compare the behaviour of the OFC model with other models, such as the BTW model or the Zhang model.

After a sufficiently long transient time, the system settles into a stationary state, where the statistical properties of the model (e.g. the probability distribution for earthquake sizes) do not depend on time. In the BTW model and in the Zhang model, the transient time is relatively brief and does not scale with system size. On average, an input energy proportional to system size is needed to reach the steady state . On the contrary, transient times in the OFC model are known to be extremely long, especially for large lattices. For example, it was claimed in Ref.  that $`4\times 10^8`$ earthquakes are not enough to reach stationarity in a system of size $`L=200`$ for $`\alpha =0.1`$. This conclusion was reached by observing the very slow convergence of the mean earthquake size to an asymptotic value during the transient time. In Ref.  though, a systematic investigation of the relaxation times for different $`\alpha `$ and $`L`$ was not attempted. A more quantitative approach to the problem was proposed in Ref.  by Middleton and Tang (MT). According to these authors a “self-organised” region develops first close to the boundaries and propagates thereafter into the bulk of the system. The distance from the boundaries of the invasion front grows with time as a power law, $`d(t)t^{\gamma (\alpha )}`$, with $`\gamma =0.23\pm 0.08,0.63\pm 0.08`$ for $`\alpha =0.07,0.15`$, respectively. The system reaches stationarity when the SOC region crosses the whole sample ($`d(t)L`$). Assuming that the power law behaviour of $`d(t)`$ holds till saturation, than the transient time of the system should scale as $`t_{tr}L^{1/\gamma (\alpha )}`$. More recently, it has been suggested in Ref.  that two distinct relaxation times exist in the system, associated respectively with the power law region and the “tail” (induced by finite size effects) of the distribution of earthquake sizes. According to this study, the former should stabilise much faster than the latter.

## III Results

Most of the studies on the OFC model at stationarity have concentrate on the probability distribution of earthquake sizes, $`P_L(s)`$, where $`L`$ is the size of the system and $`s`$ is the total number of sites that topple during an earthquake. This probability distribution does not show simple finite size scaling, at least in the range of lattice sizes accessible to simulations at present . In a recent paper , we have focused instead on the properties of earthquakes confined within a fictitious subsystem of linear size $`\lambda `$ (see fig. 1). The model is driven according to its usual dynamics but only those particular earthquakes that are entirely contained within the subsystem are counted. We have shown that if $`\lambda `$ is sufficiently smaller than $`L`$ the size distribution for this subset of earthquakes, $`P_{conf}(\lambda ,s)`$, obeys ordinary finite size scaling, i.e. $`P_{conf}(\lambda ,s)\lambda ^\beta f(s/\lambda ^D)`$, where the exponents $`\beta =3.6`$ and $`D=2`$ are universal over a range of values of $`\alpha `$.

In this work we want to address the issues briefly summarised in the previous section concerning the approach to stationarity in the OFC model. In order to be able to formulate scaling hypothesis and make use of collapse plots, we will proceed in a way similar to that of Ref. . We will consider earthquakes localised within given subsystems, in particular subsystems placed (a) at the boundaries and (b) at the centre of the system. As it is a prohibitive task to determine the time evolution of the entire distribution $`P_{conf}(\lambda ,s)`$, we will restrict ourselves to the mean earthquake size, $`<s>_{\lambda ,L}(t)`$, and, in general, to $`q`$-th moment (up to $`q=4`$) of the distribution, $`<s^q>_{\lambda ,L}(t)`$. To determine numerically these quantities we have run several simulations with different initial conditions, partitioning the time into bins of size $`\mathrm{\Delta }t`$. Let $`n`$ be the number of earthquakes occurring between time $`t\mathrm{\Delta }t/2`$ and $`t+\mathrm{\Delta }t/2`$ in a given realization of the system and let $`s_1,\mathrm{},s_n`$ be the sizes of these earthquakes. Then we define

$$<s^q>_{\lambda ,L}(t)\frac{<\left[s_1^q+\mathrm{}+s_n^q\right]_{t\mathrm{\Delta }t/2}^{t+\mathrm{\Delta }t/2}>}{<\left[n\right]_{t\mathrm{\Delta }t/2}^{t+\mathrm{\Delta }t/2}>}$$
(3)
where $`<\mathrm{}>`$ denotes an average over different realizations of the system, i.e. over different initial conditions. For each system, the parameter $`\mathrm{\Delta }t`$ has been choose small with respect to the transient time $`t_{tr}`$ but large enough to collect reasonably accurate statistics.

We consider first the case of subsystems placed adjacent to a boundary of the system, in a symmetric position with respect to the corners. In figure 2 we report $`<s>_{\lambda ,L}`$ as a function of time for $`\alpha =0.18`$ and some $`\lambda `$ and $`L`$. We observe that if the linear dimension $`\lambda `$ of the box is sufficiently smaller than the linear dimension of the system $`L`$ (approximately $`L4\lambda `$) then the curve $`<s>_{\lambda ,L}(t)`$ becomes indistinguishable for different $`L`$. This has been verified also for other values of $`\alpha `$ and $`\lambda `$. We will therefore denote with $`<s>_\lambda (t)`$ the mean earthquake size in this limit. It is already visible from figure 2 that the relaxation time of $`<s>_\lambda (t)`$ increases with $`\lambda `$. This is an indication in support of the scenario proposed by MT. Regions close to the boundaries reach stationarity sooner, signalling that an invasion front is moving toward the bulk of the system. Some of MT conclusions nonetheless will have to be modified as we will show later.

In order to describe quantitatively the invasion from the boundaries of the self-organised region we make the following simple scaling hypothesis

$$<s>_\lambda (t)=\lambda ^\eta F(t/\lambda ^z)$$
(4)
where $`\eta `$ and $`z`$ are two suitable critical exponents. In particular $`z`$ is a dynamical critical exponent that should satisfy $`z=1/\gamma `$, where $`\gamma `$ is the “invasion” exponent as defined by MT. In the limit of $`t\mathrm{}`$ the scaling function $`F(x)`$ saturates to a constant, implying that the exponent $`\eta `$ is related to the finite size exponents of the probability distribution $`P_{conf}(\lambda ,s)`$ by the relation $`\eta =2D\beta 0.4`$. In figure 3 we report collapse plots of the form (4) for various values of $`\alpha `$. We observe that a reasonably good collapse could be obtained for all the $`\alpha `$ if we choose the universal exponent $`z=1.3\pm 0.1`$. We are therefore led to conclude that a universal exponent $`z`$ exist contrary to the claims by MT.

We consider next subsystems of different sizes placed at the centre of the system. We observe that the relaxation time does not depend on the size of the subsystems (see fig. 4a). This confirm that the self-organisation mechanism develops from the boundaries and that the system enters stationarity when the self-organised region span the whole system. Only when the bulk of the system is reached by the self-organised region is stationarity settled so that concentric subsystem of different sizes will inevitably reach stationarity at the same time. We have also tested the relaxation times for higher moments of the avalanche probability distribution to see whether different parts of the distribution (e.g. power law part and the “tail”) have different relaxation times as suggested in Ref. . In our investigation though we have not observed any significant difference in the relaxation times associated with different moments. We report as an example in fig. 4b the behaviour of the first, second and fourth moment in a particular case.

The scaling equation (4) suggests that the transient time in the OFC model scales with system size as $`t_{tr}L^z`$. One way to test this is by comparing the time behaviour of the average earthquake size in a central subsystem of size $`\lambda `$ for different system sizes $`L`$. Indeed the asymptotic value $`<s>_{\lambda ,L}(t\mathrm{})`$ should not depend on $`L`$. We report as example the case for $`\alpha =0.18`$ in fig. 5, where we have rescaled the time by a factor $`L^z`$. The curves show some noisy behaviour, due to the difficulties in collecting good statistics (relatively few earthquakes occur in the bulk of the system as compared to the boundaries). Nonetheless the value deduced for the exponent, $`z1.3`$, is consistent with the determination made through the analysis of the earthquakes occurring at the boundaries. We have obtained similar results also for other values of $`\alpha `$. In addition, besides the average earthquake size, we have considered also the time behaviour of other quantities such as the roughness of the energy landscape (in analogy to surface growth problems) and the number of earthquakes per unit time. All these different quantities on average reach stationary at the same time.

The algebraic divergence of the relaxation time with system size reflects the presence of long-range spatial correlations in the stationary state. Indeed if correlations were only short range, than one would expect that the transient time would not scale with system size. This is for example the case for the BTW model in $`d`$ dimensions, where the height-height correlations are algebraic but decays as fast as $`r^{2d}`$ ($`r`$ being the distance between two sites) . We have measured for various $`L`$ and $`\alpha `$ the probability distribution of the spatially averaged force in the system

$$M=\frac{1}{L^2}\underset{i,j=1}{\overset{L}{}}F_{i,j}$$
(5)
In a system with sufficiently short range correlations, this probability distribution would tend, in the limit $`L\mathrm{}`$, to a gaussian distribution around the mean due to the central limit theorem (this is indeed what results in the BTW model). Let $`<M>`$ and $`\sigma _M`$ be respectively the average and the standard deviation of the distribution. In figure 6 we have plotted $`\mathrm{log}(\sigma _MP(M))`$ versus $`(M<M>)/\sigma _M`$ for various $`L`$ and $`\alpha `$. Using these coordinates a gaussian function would result in an inverse parabola. For each $`\alpha `$ value the data of figure 6 collapse on a single function, which is clearly not gaussian (deviations from gaussianity are more pronounced for increasing $`\alpha `$ values). This indicates that the central limit theorem does not hold in this case, not even for large $`L`$, suggesting that long range algebraic correlations are present and therefore the sum (5) can not be decomposed into a sum of independent terms. This observation is in agreement with the results reported in Ref.  where the presence of long-range spatial correlations were deduced from the behaviour of a suitably defined susceptibility, $`\chi (L\sigma _M)^2`$. It was claimed that $`\chi `$ diverges as $`L^2`$ and correspondingly that $`\sigma _M`$ is, to leading order, independent on $`L`$ (if $`M`$ was a sum of uncorrelated variables, $`\sigma _M`$ would decrease as $`1/L`$). In our investigation we have found that $`\sigma _M`$ slightly increases with $`L`$ but asymptotically tends to a constant value, in accordance with Ref.  ($`M`$ is a bounded variable so $`\sigma _M`$ can not grow indefinitely).

## IV Conclusions

In conclusion, in this paper we have examined the process of self-organisation in the OFC model. By considering earthquakes confined within a given subsystem we have been able to clarify some of the issues related to this problem. In accordance with Middleton and Tang  we have found that SOC develops first close to the boundaries and subsequently invades the interior of the system. The invasion process is controlled by a dynamical critical exponent, $`z1.3`$, which, contrary to previous claims, is universal over a range of values of the dissipation level of the local dynamics. We have shown moreover that the transient time in the system scales with system size as $`t_{tr}L^z`$. This is a peculiarity of the OFC model as other “sandpile-like” models (e.g. the BTW and the Zhang models) do not display any scaling in the transient time. We have associated this feature with the presence of long-range spatial correlations in the stationary state.

Our findings are in general agreement with recent works on the OFC model . Indeed we have provided complementary evidences (not based on the probability distribution for earthquake sizes) that the model is critical even in a non conservative regime. Moreover it confirms that there is universality in the system and that finite-size scaling can be recovered by considering subsystems whose linear extent is sufficiently small.

Finally, it is interesting to remark that the probability distribution for the spatially averaged force in the system is somewhat reminiscent of a probability distribution observed in a confined turbulent flow experiment  (BHP). As a term of comparison we have reported in fig. 6 the BHP functional form over-imposed to the curve for $`\alpha =0.21`$. Attempts to link SOC systems to turbulent phenomena have long been suggested, but only recently this has been put on a firmer basis .

We thank M. Paczuski for helpful conversations and H.J. Jensen for a critical reading of the manuscript. This work was supported by the EPSRC (UK), Grant No. GR/M10823/01 and Grant No. GR/R37357/01.


format: text

text:
# Chandra Observation of LMXBs in the elliptical galaxy M84

## 1.  Introduction

Observation of X-ray sources in the Milky Way galaxy began with the launch of the first X-ray detector (Giacconi et al. 1962). The study of accreting systems at all wavelengths is now a well-developed field (for a recent review see Tanaka & Shibazaki 1996; theoretical progress is presented by Iben et al. 1995; Kalogera & Webbink 1998). With the advent of high angular resolution X-ray observations, first with Einstein and ROSAT and now with Chandra and XMM-Newton, it has become feasible to study the X-ray populations in other galaxies.

This Paper presents a detailed study of Chandra observation of the point sources in M84, an elliptical galaxy in the core of the Virgo cluster (study of the diffuse X-ray emission of M84 is presented in Finoguenov & Jones 2001). The particular importance of population studies in early-type galaxies comes from the small number of possible avenues leading to the X-ray source production. As most ellipticals are old stellar systems, only accretion of matter in low mass X-ray binaries (LMXB) can provide enough power to explain X-ray luminosities that exceed $`10^{37}`$ ergs/s. The question of how low is the mass of the secondary star is closely connected to studies of the stellar population (e.g. Worthey 1994) and recent merger activity (e.g. Schweizer & Seizer 1992). The latter can be responsible for differences in the X-ray population in the centers of galaxies, compared to that in the halos. The dependence of LMXB activity on star-formation rate is considered in White & Ghosh (1998).

When one star in a binary system becomes a compact object (primary), the accretion process can be triggered through the Roche-lobe overflow of material from the evolved low-mass companion star (secondary). Another mechanism for triggering accretion at epochs long after the cessation of star-formation, is by shrinking the binary star separation through magnetic braking (outward transport of the angular momentum in magnetic stellar winds, Verbunt & Zwaan 1981) or gravitational radiation losses (Paczyński 1967) leading to accretion even from a main sequence star. Accretion onto a white dwarf is identified through a super-soft X-ray spectrum (e.g. Kahabka & van der Heuvel 1997). However, our observation is not sensitive to this mechanism, since detecting white dwarfs in our softest energy band ($`0.31`$ keV) would imply a luminosity for the source (calculated assuming a spectral temperature of 35 eV), far exceeding the corresponding Eddington limit.

The luminosity function of X-ray sources in M84 can be expressed as a convolution of the accreting binary mass-function and X-ray luminosity function for equal mass primaries. Neutron stars, the most numerous class of accreting binary primaries, have very similar masses, 1.4 M. Thus, the source luminosity function below the Eddington luminosity for a neutron star ($`1.8\times 10^{38}`$ ergs/s) is dominated by the X-ray luminosity function of neutron-star LMXB (NS LMXB). Their luminosities are proportional to the accretion rate determined by the orbital period and mass loss (Wu 2001). Most of the known black-holes in binary systems have a mass of 6 M (Tanaka & Shibazaki 1996), however recent studies indicate a more complicated black hole mass function (Bailyn et al. 1998; Fryer & Kalogera 2001).

In this Paper we study the spatial distribution of LMXB in M84 (§3) and derive the integral $`log(N)log(S)`$ distribution (§4). We study source hardness ratios and conclude with a discussion on the emission mechanism for LMXB (§5). Interpretation of the M84 LMXB luminosity function is given in §6.

## 2. Analysis of the ACIS observation

In the ACIS-S 28.7 ksec observation of M84, point source detection was done in the 0.4–5.0 keV energy band, where the source signal-to-noise is highest. Throughout sections §3–§5, luminosities are cited for the 0.4–10 keV band, using the countrate-to-luminosity conversion for the best-fit power law photon index of 1.4, derived for the combined spectra of all point sources within two effective radii (Finoguenov & Jones 2001). Given the duration of the observation, and assuming a 17 Mpc distance to M84, one ACIS-S count corresponds to $`10^{37}`$ ergs/s. Estimates of the detection rate for the background sources were made using a power law photon index $`\mathrm{\Gamma }=1.7`$ and the $`log(N)log(S)`$ determination by Giacconi et al. (2001).

Source detection was carried out using the matched filter (wavelet) technique (Vikhlinin et al. 1995). We performed simulations to estimate the detection efficiency of this method on our field, as described below (§4). In addition we estimate the effects of variations in the PSF on the resulting signal-to-noise and thus source detection efficiency. Variations in the PSF cause the apparent source size to increase with increasing off-axis angle. Thus, we run the detection procedure twice, first for the central $`10^{}`$ radius region, where we use a full resolution ACIS pixel image ($`0.492^{\prime \prime }`$ on a side) and second employing a 4 by 4 pixel binning for the full $`15.6^{}\times 11.1^{}`$ field, shown in Fig.1, which includes S3 and portions of S2, I2 and I3 chips. At large radii, a decrease in the detection efficiency arises due to the larger detection region for the source counts and therefore for the background. We merged the source lists from the two methods, thus removing duplications. To determine the source count rate, we calculate exposure maps, accounting for the spatial non-uniformity of the CCD quantum efficiency. We take vignetting (weighted by the mean spectrum of sources, $`\mathrm{\Gamma }=1.4`$) and dithering into account, using CIAO contributed software.

## 3. Spatial distribution of X-ray binaries in M84

Before we proceed to describe our results, we note that for M84, we can study the X-ray population at a few effective radii, while strong diffuse emission prevents us from detecting faint sources very close to the galaxy center.

The Chandra ACIS image of the M84 field is shown in Fig.1. We show X-ray images of the full field and the center of M84 with detected sources marked with circles. The enhanced background (in the upper left of the top left image) is from a neighboring CCD chip (S4), which we omitted from analysis.

Of the 106 sources shown in Fig.1 we exclude 3 sources associated with foreground objects from all subsequent analysis. We identified these as corresponding to an optical counterpart of size or luminosity exceeding that expected for globular clusters at the M84 distance. Spatial analysis of the source density does not reveal any peculiarities at the 5$`\sigma `$ level. At the 3$`\sigma `$ significance level (Fig.3), a clustering of sources on sub-arcminute angular scale is seen, as well as the effect of non-uniform source detection in regions with strong diffuse emission, corresponding to the $``$-shape of the hot gas in M84, seen in the upper panels of Fig.1. Strong clustering on sub-arcminute angular scales is an attribute of the sources that constitute the Cosmic X-ray Background (e.g. Vikhlinin & Forman 1995).

The profile of the source number density, centered on M84, is presented in Fig.3. While the outer part is described well by the stellar light distribution, taken as a de Vaucouleurs (1/4) law with $`r_e=1^{}`$ (where $`r_e`$ is the effective radius, within which the galaxy contains half of its light), the center reveals a deficit of sources, due to the high detection threshold, caused by the bright diffuse X-ray emission. In the next section, we develop a method to correct for this effect, with the results also shown in Fig.3.

With the normalization shown in Fig.3, and assuming a distance to M84 of 17 Mpc (the corresponding luminosity in the B band is $`L_B=4.5\times 10^{10}L_{}`$ within $`r_e`$), the LMXB rate per $`L_B/10^{10}L_{}`$ is $`23\pm 2`$ for X-ray luminosities $`>3\times 10^{37}`$ ergs/s. The fact that the X-ray binary distribution is similar to that of the optical light, strongly suggests a similar formation history for single and binary stars. Latter mixing should not play a strong role, since metallicity gradients are detected in this galaxy (Kobayashi & Arimoto 1999), and otherwise would be erased.

## 4. $`log(N)log(S)`$

In deriving the efficiency for source detection, we simulated the effect of the diffuse X-ray emission on the resulting signal-to-noise requirements for our source detection at a $`5\sigma `$ statistical threshold. We also took into account the spatial distribution of sources, which are concentrated toward the center, where the ’background’ is higher. This further reduces the detection sensitivity. These effects are important for source detection close to the center of M84. For sources outside the effective radius of M84, changes in the sensitivity (due primarily to quantum efficiency differences between CCD chips) and the reduced signal-to-noise, due to PSF degradation at large off-axis angles, play important roles in the detection. We estimate the resulting detection efficiency in Fig.4. The effects of strong diffuse emission are important for sources with fewer than 15 counts, while changes in the PSF affect sources with 6 counts or fewer. Because of the Poisson noise, the shape of the source luminosity distribution will be further modified (Poisson migration) and will be different for different exposure times. While a proper reconstruction requires a response matrix with the luminosity resolution defined by a Poisson process, in the present analysis, we simply correct for the amplitude of the effect, by convolving the detection efficiency with the Poisson probability distribution and show the resulting detection efficiency (which is equivalent to the completeness) in Fig.4. In the following derivation of the $`log(N)log(S)`$, we will correct for this incompleteness in the survey. Finally, changes in the sensitivity (effective area and quantum efficiency) are most important for flux correction and affect the resulting slope at high fluxes.

The resulting $`log(N)log(S)`$ is shown in Fig.4. In further analysis we accounted for the contribution from extragalactic sources (Giacconi et al. 2001). There is a break at the level of 19 sources (with an expected CXB contribution of 6 sources) at a corresponding luminosity of $`L_b=4.6_{0.6}^{+1.2}\times 10^{38}`$ ergs/s (errors quoted are 68% confidence limits). The integral slope at lower fluxes is $`0.87\pm 0.07`$, while at higher fluxes, it is $`1.8\pm 0.5`$. The differential slopes at lower and higher fluxes are $`1.79\pm 0.18`$ and $`2.7_{1.1}^{+0.6}`$, correspondingly. Our results for the luminosity of the break and the slope at high fluxes are in remarkable agreement with findings by Sarazin et al. (2000) for another elliptical galaxy, NGC4697, ($`L_b=3.2_{0.8}^{+2.0}\times 10^{38}`$ ergs/s and differential slopes at lower and high fluxes of $`1.29_{0.36}^{+0.49}`$ and $`2.76_{2.0}^{+0.8}`$, correspondingly). However, as we will argue below, luminosities of the break in both galaxies are derived for the averaged spectrum, which is not typical for NS LMXB emitting close to the Eddington.

Flattening of the source luminosity function below the luminosity of $`5\times 10^{37}`$ ergs/s, seen in NGC4697 and M31 also may be present in M84, but since we have only one point there, derived from our detection of sources with only 3 and 4 counts, it would certainly be an over-interpretation to discuss this in more detail.

## 5. Hardness ratios

In the absence of a variation in absorption across M84, plus the lack of neutral gas in ellipticals, the source hardness ratios are particularly straightforward to interpret.

In Fig.2 we present the hardness ratios for the 32 sources with the best statistics, determined using the source counts in the Hard and Soft bands as $`(HS)/(H+S)`$, using three independent energy bands, $`0.31`$, $`12`$ and $`210`$ keV. Given the precision of the hardness ratio estimates, a combination of first vs second (HR21), and second vs third (HR32) reveals the largest spread in the points and therefore provides better tests for the predictions from various spectral models.

We subdivide the source sample into two parts, within or outside two effective radii, as we expect almost no background sources to be detected in a small region close to the center of the galaxy. The larger number of sources available at larger radii, allows us to better sample the X-ray colors. The $`log(N)log(S)`$ of M84 sources and of cosmic X-ray background sources are quite similar in the observed flux range. In the whole field used for source detection, CXB sources constitute almost 1/3 of the detected sources, i.e. in Fig.2, 11 background sources should appear, including 6 of high statistical significance (solid circles). Given the energy band of Chandra and using the results of point source hardness ratios and $`log(N)log(S)`$ analysis from the deep XMM pointing of the Lockman hole (Hasinger et al. 2001), we conclude that 1/3 of the background sources should exhibit hard colors (within the power law – nH grid shown in Fig.2, note that a variation of nH is intrinsic to AGN), while the other 2/3’s should be within the grey-shaded area on the right panel of Fig.2 of ROSAT-type background sources (type 1 AGN). In fact a fraction of the ROSAT-type background sources should be larger, due to differences in the sensitivity between Chandra and XMM at high energies.

We find two concentrations of hardness ratios. The first (black circles in Fig.2) corresponds to that expected for a multicolor disk black-body model (Makishima et al. 1986; diskbb model in XSPEC) with a central black-body temperature near 0.5 keV. The second (gray circles in the shaded area on the right panel of Fig.2) corresponds to power law indices between 1 and 2. The empirical relation between the luminosity and the spectral shape, derived for the Milky Way sources, implies that luminous high state LMXB’s should be dominated by a thermal component, while for sources in the low state, the spectrum is harder. This is in a good agreement with our detection of a number of sources with thermal disk emission, which includes most of the sources with flux levels near the break in the luminosity function. Thus, although the disk black-body model cannot account for all the observed X-ray colors, the large concentration of M84 sources about the prediction of this model allows us to distinguish them from CXB sources. Among the remaining sources there are four that follow the prediction of the black-body spectrum (the furthest to the left dotted line in the comptonization grid in the left panel of Fig.2). Their luminosities exceed the Eddington limit for a neutron star. Other sources exhibit both a soft-excess and a hard tail (these sources fall in the lower right of the diagram). The hard tail in M84 sources could be produced by comptonization on the hot plasma ($`kT10`$ keV) with typical $`\tau =13`$ and a soft input spectrum $`kT<0.1`$ keV (Sunyaev & Titarchuk 1985; Titarchuk 1994; comptt model in XSPEC).

Our conclusions regarding the origin of X-ray colors in M84 LMXBs are in agreement with a much more developed scheme for X-ray binaries in the Galaxy (e.g. Done 2002). The next step would be then to employ the luminosity-spectrum relation, established for the Milky Way sources to separate the NSs from BHs in M84 at low luminosities. Given that the luminosities of M84 sources with well measured colors exceed $`10^{38}`$ ergs/s, the soft component should still be dominant in LMXBs with a black-hole primary, yet their Milky Way prototypes exhibit a large spread in colors. On the other hand, the majority of LMXBs with NS-primary should be Z-sources (like Cyg X-2), that are characterized by a modest variation in colors. Following these guidelines, X-ray black holes comprise up to 1/3-1/2 of the X-ray sources in M84. Chandra observations of the Sombrero galaxy (Delain et al. 2002), indicates a somewhat smaller, yet similar fraction of accreting black holes (1/4). While small-number statistics prevents us from a firm conclusion, the observed differences are in line with the suggestion of a gradual built-up of black-holes in the 1.5-3 M range (Prokhorov & Postnov 2001).

## 6. Discussion

### 6.1. The break: Eddington limit on isotropic radiation from accretion of matter on the neutron star.

We note in §5, that sources near the break of the luminosity function have similar spectra. Their spectra are softer than the average one. Therefore, using a softer spectrum we find the break in the luminosity function occurs at $`L_b=2.4_{0.3}^{+0.6}\times 10^{38}`$ ergs/s (a factor 0.52 lower than the estimate using the averaged spectrum). The Eddington limit also should be corrected for the $`\mu _p=<A>/<Z>`$, which for the typical solar metallicity of stars in M84 is 1.17. Thus the corrected value for the isotropic emission from accretion of matter of solar elemental composition onto a neutron star is $`2.1\times 10^{38}`$ ergs/s, which agrees well with our measurements of the break.

The existence of a limiting luminosity arising from accretion onto a neutron star has been debated by several authors. Paczyński & Wiita (1980) proposed that at supercritical accretion rates, the resulting luminosity could be much higher than the Eddington luminosity. Thus, our result of a very close match between the break luminosity and the Eddington luminosity is important in modeling the structure of the inner accretion zone. Inogamov & Sunyaev (1999) find that at accretion rates close to critical, emission from the surface of the neutron star becomes softer, comparable to the emission from the disk. Compared to emission from the disk, spectrally harder emission from the surface of the neutron star is expected at lower accretion rates and together with black holes is probably responsible for the harder colors of the combined spectrum.

A prediction of a cut-off in the luminosity function is not unique to the Eddington limit. For example, on a basis of LMXB evolution, Wu (2001) predicts a gradual decrease in the number of accreting systems with high mass transfer rates with galactic age, with a luminosity function exhibiting an exponential cut-off. However, a close correspondence of the break in the luminosity function to the Eddington luminosity, also seen in star-forming galaxies (Körding, Falcke, Markoff 2002), is a strong argument in favor of the relation of the break to the limiting luminosity.

### 6.2. The high end of the luminosity function: the mass function of accreting black holes

Chandra observation of M84 reveal 5 to 10 sources (the number is sensitive to the behavior of the luminosity function just above the break), whose luminosity exceeds the Eddington limit for accretion onto a neutron star. When BH LMXBs (LMXBs with a black hole primary) are used to explain the high-luminosity end of the source population in M84, a problem arises in explaining why such a steep dependence of source number on luminosity becomes flatter below $`2\times 10^{38}`$ ergs/s. One can seek an explanation in terms of the specific BH LMXB luminosity function111We define the specific luminosity function, as a distribution of luminosities (resulting from a distribution of the accretion rates) for a given mass black hole., but such a coincidence with the Eddington limit is unlikely. On the other hand, it is natural to attribute this to the black-hole mass-function. In this case, the differential slope of the black hole luminosity function should be shallower than -1.79, or no break would be seen, while a steeper slope of the luminosity function above the break is a result of convolution of the black hole mass and luminosity functions. Using our measurements, we can constrain the differential slope $`\alpha `$ of the mass function for accreting black holes in M84 to be $`2.7<\alpha <0.9`$, where the lower limit comes from an assumption that black holes shine at their Eddington limit and the upper limit corresponds to similar luminosity functions of LMXBs with black-hole and neutron star primaries. The statistical significance of the slopes quoted above is the same as for the slope at high luminosities ($`\genfrac{}{}{0pt}{}{+0.8}{2.0}`$).

Besides accreting black holes, explanations for the appearance of super-Eddington sources, based on an accreting neutron star could be: beaming of the radiation, flaring and super-Eddington accretion rates. Each of these scenarios has difficulties explaining the observed luminosity function of M84. Observational examples of super-Eddington accretion rates, such as for SS433, demonstrate that the luminosity does not exceed the Eddington limit by more than a factor of two. Occurrence of a break in the observed luminosity function at the Eddington limit for isotropic radiation would be hard to reproduce in the beaming scenarios. As the existence of stellar mass black-holes is not disputed, our conclusion on the effect of black-hole mass function provides most likely explanation for the observed behavior of the M84 source luminosity function at high fluxes.

The brightest candidate galactic source in M84 has a luminosity of $`1.6\times 10^{39}`$ ergs/s. In the inner 2 region of M84, where most of the galactic sources should be located and where the CXB contamination is low, the brightest source has a luminosity of $`0.9\times 10^{39}`$ ergs/s. This is contrary to the very luminous sources found in the Antennae galaxies and other star-forming galaxies, whose spectral characteristics also are different (Fabbiano, Zezas, Murray 2001; Makishima et al. 2000). Interpretation of the Antennae sources could be either via beaming (King et al. 2001; Körding, Falcke, Markoff 2002; Zezas & Fabbiano 2002) or exceeding the Eddington limits (Begelman 2002).

### 6.3. Low end of the luminosity function: sampling the initial distribution of binaries with NS primaries

To relate the observed luminosity function to the distribution of objects with different physical parameters, we need to make an assumption about the origin of the accretion process. Two major scenarios are possible, accretion from a main sequence star or accretion caused by nuclear evolution of the secondary.

The first scenario dominates for LMXBs in the Milky Way and can accommodate both the high accretion rates and short periods of the LMXB systems. To apply this model to the observed luminosity function of LMXBs in M84, we need to consider systems with relatively high accretion rates (typically in excess of $`10^9`$ M yr-1). At these accretion rates, assuming they are persistent, observed LMXBs must have formed recently (less than a Gyr ago). The old age of the stellar population (12 Gyr) implies that the fraction of active systems is 1/100–1/10 of the integral LMXB production (a factor of $`\mathrm{\Delta }T/T`$). Yet, as the time since star formation increases, different mechanisms come into play. However, for accretion from a main sequence star, magnetic braking (Verbunt & Zwaan 1981) can change the behavior of period over time. For example, for a similar mass range, which we take as $`0.31`$ M for the secondary, and using the prescription for the magnetic braking from Wu (2001), at 12 Gyr after star-formation, binaries with longer initial periods (0.3-1.3 days) are important, compared to those with periods shorter than 0.3 days that dominate the accretion process only 0.1 Gyr after star-formation. Therefore, the LMXB population in galaxies of different ages essentially samples different parts of the initial period distribution. The effect should be seen in comparison studies between galaxies and will provide clues to binary formation (Kalogera & Webbink 1998).

For LMXBs with an evolved secondary, the luminosity and period of the system are not independent. Assuming a 20% efficiency for accretion, the luminosity of LMXBs is given by $`0.4\times 10^{37}P_d^{0.93}M_2^{1.47}`$ ergs/s for accretion from an evolved companion, where $`P_d`$ and $`M_2`$ are the period of the binary in days and mass of the secondary in M (Webbink et al. 1983). Given the typical solar mass for the secondary, the luminosity of $`10^{38}`$ ergs/s requires periods on the order of ten days. Most LMXBs in the Milky Way have short periods, but the Milky Way star-formation also could be recent (note that the observational appearance of LMXBs in the Milky Way bulge is generally attributed to high kick velocities, e.g. White & Ghosh 1998). While studies of periodicities of LMXBs in ellipticals are the most direct way to disentangle systems with evolved secondaries from those with main sequence secondaries, the required collecting area for such observations will only be available with new generations of X-ray telescopes (e.g. XEUS). Element composition of the accreting material from an evolved secondary is dominated by helium. This increases the Eddington limit by a factor of two. Therefore, observing a lower value for the Eddington limit, favors accretion from main sequence stars in LMXBs.

Finally, we would like to comment on the production of LMXBs in globular clusters. Stellar trapping by a compact object in globular clusters can lead to LMXB formation (Kuranov, Postnov, Prokhorov 2001). Although the number of X-ray sources associated with globulars exceeds the average expected for their light, the number is still a small fraction of the total number of sources ($`10`$%) in some ellipticals (Sarazin et al. 2000; Kraft et al. 2001). Therefore, if globulars are a primary site for LMXB production in early-type galaxies, to explain the LMXBs outside globular clusters either globular clusters have been disrupted or a “kick” expelled the LMXB (as during the collapse of the white dwarf into the neutron star). Indirectly the importance of globulars can be estimated by comparative studies (White 2001). In addition, the outskirts of brightest cluster galaxies, for example NGC1399 and M87, are characterized by an astonishingly high frequency of globulars, which explains a strong association of globulars with X-ray sources in the outskirts of NGC1399 found by Angelini, Loewenstein, Mushotzky (2001).

## Acknowledgments

This work was supported by NASA grants GO0-1045X and AG5-3064 and the Smithsonian Institution. AF has benefited from discussions with Andrei Beloborodov, Eugene Churazov, Guenther Hasinger, Vicky Kalogera, Wolfgang Pietsch, Konstantin Postnov, Juri Poutanen, Sergey Sazonov, Rodrigo Supper, Yasuo Tanaka and Andreas Zezas at various stages of the presented research. AF acknowledges receiving the Max-Plank-Gesellschaft Fellowship. Authors thank the referee for useful comments on the manuscript.


text:
# Stress transmission through three-dimensional ordered granular arrays

## I Introduction

In a static granular bead pack all of the forces on a particle, due to the contacts from its neighbors and due to gravity, must be in perfect balance; there can be no net force on any particle. Moreover, in the limit of hard particles, no particle distortion is allowed leading several investigators to question whether elasto-plasticity theory can be used in this hard-sphere limit to describe the mechanical properties of an amorphous granular material . Others suggest that elastic theory is capable of producing the observed experimental results . The presence of disorder makes an exact calculation of the contact forces in such a material impossible. On the other hand, it should be possible to calculate the forces exactly in a perfect crystal . However, because macroscopic grains do not spontaneously crystallize, it has been difficult to create large defect-free granular crystals  on which to perform an experiment to see how the forces are transmitted. In addition, there will always be small imperfections or asperities in the particles themselves and slight variations in the diameters of the grains. Although we cannot completely avoid disorder of this latter kind, we have found an efficient method of constructing essentially perfect crystals with few, if any, defects in the particle positions . With such crystals, one would hope that the presence of small disorder from the variation of the individual particle shapes and sizes could be treated as a perturbation about the perfect crystalline response.

Previously, there have been experimental studies of the statistical properties of the local contact forces in granular packs. In particular, the distribution of normal forces, $`P(F)`$, between the particles and the container walls was measured in a variety of situations . In most of those studies, a piston pressed down on the entire top surface of the granular material. The surprising result was that the entire form of $`P(F)`$ (a large value near $`F=0`$, a small peak \[or plateau\] with a maximum near the average force and an exponential tail at large forces) was robust and did not depend on whether the granular medium was amorphous or crystalline or whether the particles were smooth or rough . The question is thus raised as to what aspects of the force propagation in a granular material are affected by preparation and the underlying crystal or amorphous structure.

Several experiments have studied the response (i.e., the Green’s function) to a localized external force at the top surface of a granular medium. The different studies found somewhat contradictory results. DaSilva and Rajchenbach  studied a system consisting of a special two-dimensional packing of photoelastic rectangular bricks. They measured the contact forces between particles as a function of position and found that the contact forces had a maximum centered below the position of the applied perturbation and that on ether side the forces decreased in magnitude. The width of this peak in the forces grew as the square root of the depth below the surface. This result is consistent with models predicting a diffusive propagation of the forces such as the q-model proposed by Coppersmith et al. . Other studies have focused on the response to a localized external force in amorphous packings of spheres. Reydellet and Clément found that the response at the bottom of a three-dimensional pack had a central maximum of approximately Lorenztian shape with a width that grew lineraly with the depth of the pack . Similar behavior was found in two-dimensional samples studied by Geng et al. . This linear increase of the width with depth is consistent with elastic theories . The work by Geng et al. also suggests a strong effect due the spatial ordering of the particles resulting in a central minimum of the response forces in ordered packings. Thus there is a need to study the propagation of forces from a localized perturbation in a nearly perfectly ordered crystalline array of particles in order to investigate how the underlying order (or disorder) of the lattice influences the force propagation.

Bouchaud et al. have proposed a model in which forces propagate in straight lines until encountering a defect where they are split, and they have shown that this model leads to a central force minimum for small depths and to elastic-like behavior for large pack depths . The straight-line propagation of forces (resulting from the hyperbolic model of Cates et al. ) and the splitting of forces can be viewed in a geometrical framework based upon the contact orientations of neighboring beads. From this viewpoint crystalline packings allow for the control of these contact orientations and the resulting splitting. In particular, the hexagonal-close-packed (hcp) crystal structure is such that forces are required to split at nearly every layer as they travel downward through the pack. By looking at the response of an hcp crystal to a localized force it is possible to explore the consequences of multiple splittings in comparison with the predictions of the forcechain splitting model.

In this paper we present an experimental investigation of the response to a localized external force in three dimensional nearly perfect hcp and face-centered-cubic (fcc) crystals. We measure the resulting contact forces at both the top and bottom surfaces of the crystal, and investigate how the pattern of forces depends on the crystal structure and the depth of the pack. We determine that in the regime of highly ordered granular packs, the crystal structure plays a dominant role in determining the average pattern of response forces at the bottom surface. These patterns can be viewed as a consequence of force balance in the specific geometry of the crystal. We also consider the effect of small amounts of disorder in these crystals and show that they modify the pattern of forces on the bottom surface in a nontrivial manner and scatter forces back upwards through the pack to the top surface.

## II Experimental Methods

We created large nearly perfect hcp and fcc granular crystals of $`3.06\pm 0.04`$ mm diameter soda lime glass spheres. These crystals of approximately $`53`$ bead diameters on a side were contained within an acrylic cylinder with a triangular cross-section designed to be commensurate with the hcp crystal structure. The crystals were created as individual layers of triangular order in the manner described in reference . In this way the position of each individual bead was viewed and corrected if neccessary as the crystal was constructed, so that the resulting crystal would have very few defects (estimated to be less than ten beads significantly out of place for a crystal of approximately $`50000`$ beads). Furthermore, the stacking order of the planes was controlled in order to achieve three dimensional fcc and hcp crystals.

After building such a crystal we applied an external force to a localized region at the center of the top of the crystal either by using a slowly applied steady force of approximately $`100`$ pounds or by using a quick impulse force. When using the slowly applied force it was necessary to apply the force over a region approximately six beads in diameter. Attempts to apply the force to a smaller region resulted in significant bead movement. When using the quick impulse force it was possible to apply the force to a smaller area of approximately two beads in diameter. For the purpose of comparison experiments were also done with the quick impulse force applied to the large six bead diameter region.

We used a carbon paper technique to measure the forces in response to the external forces . A piece of carbon paper and a piece of white paper  were placed on the top and bottom of the crystal. As the force was applied individual beads pressed into the carbon paper and left marks on the white paper. The size and darkness of the resulting marks was related to the normal force on the corresponding beads. When using the slowly applied force, we calibrated the force versus darkness response as explained in . The positions of the carbon marks were matched to a triangular lattice and the forces for each lattice point were averaged over approximately ten experimental runs. It is important to note that individual runs vary significantly, but that the average over many runs remains highly reproducible. In order to improve statistics, the locally measured forces were also averaged over symmetries of the crystal (one reflection and two rotations).

Note that the analysis technique must be carried out carefully so as not to introduce error from the averaging over many files, and crystal symmetries. If the centers of the images are not determined accurately, then lattice positions from file to file will not match up exactly and will result in a small degree of broadening of any features. We estimate this broadening to be significantly less than one bead diameter for each of the regions of large force.

When using the quick impulse, we calibrate the measured forces by comparing them with those obtained with an identical set-up using the slowly applied force described above. We compared the patterns of carbon marks left by the two methods and found no qualitative difference. Figure 1(a) shows the pattern of calibrated forces at the bottom of a $`19`$ layer fcc crystal in response to a slowly applied force at the top surface over a region six beads in diameter. The calibration maps points of zero intensity to a small but nonzero force which results in a nonzero background. Figure 1(b) shows the pattern of intensities of carbon marks at the bottom of a similar crystal in response to a quick impulse force applied to the same area. The color scale of intensity has been shifted by a factor of ten and a minimum value used to match the previous pattern of forces. Using this shift in color scale intensity we map the two different measurements onto one another and obtain a calibration for the quick impulse.

## III Results

In amorphous packings it is difficult to apply a localized force large enough to leave carbon marks without also having significant movement of the grains. As such we were only able to explore the response of amorphous packings to the quick impulse force applied to a large six bead diameter region. In this case we observed the pattern of intensities at the bottom of three-dimensional amorphous packs of various depths and found a broad central maximum. Our data on the width of this central maximum are consistent with a square-root growth with depth as seen in a two-dimensional packing of bricks by DaSilva and Rajchenbach , but are also within error of a linear growth as seen by Reydellet and Clément  in three-dimensional packings.

### A FCC crystals

When using highly ordered packings bead movement is much less of a problem, and it is possible to use a slowly applied force to a large area or a quick impulse to a smaller area. For fcc crystals we found a characteristic response pattern at the bottom surface which is qualitatively independent of type of force applied. Figure 2 shows the patterns of intensities at the bottom of fcc crystals, $`9`$ layers and $`19`$ layers in depth, in response to a quick impulse force applied to an area two beads in diameter at the center of the top surface of the crystal averaged over approximately ten experimental runs and averaged over crystal symmetries. We find three regions of large force, arranged in a triangular pattern. The size of this triangle increases with the depth of the crystal. The size of the three regions of large force are comparable to the size of the region over which the force was applied. In addition to the three regions of large force, there appear fainter but still significant regions of force along the lines connecting these three regions.

Figure 3(a,b) shows cross sections of the intensity versus position across the triangular pattern at the bottom surface of $`9`$ and $`19`$ layer fcc crystals. As depicted in the inset of figure 3(a) the cross-section proceeds across the edge of the triangular pattern, corresponding to a small peak, through the center of the pattern (position $`=0`$) with small intensities, and across the region of large force, corresponding to a large peak. The position of the large peaks as a function of the depth of the crystal are indicated in figure 3(c), with the error bars corresponding to the full width at half maximum of the peaks, which can in some cases be seen to be significantly skewed. Additional data points are included in part c corresponding to experimental runs with a quick impulse applied to a large area and to runs with a slowly applied force over a large area. The sloped line indicates an angle of $`35.3`$ degrees which corresponds to the angle down the edge of a tetrahedron. The regions of large force are found to lie just inside of this angle as explained in section IV. The inset of figure 3(c) shows the position of the small peaks corresponding to the edge of the triangular patterns as a function of the depth of the crystal packing. The solid line shows an angle of $`19.5`$ degrees with respect to the vertical, which is corresponds to the angle down the face of a tetrahedron.

### B HCP crystals

In contrast to the triangular patterns found for fcc crystals, hcp crystals show rings of large force at the bottom surface in response to a localized force at the top surface. Figure 4 shows the intensity patterns at the bottom surfaces of hcp crystals $`9`$ and $`25`$ layers in depth in response to a quick impulse force applied to an area two beads in diameter at the top surface averaged over approximately ten experimental runs and averaged over crystal symmetries. As the depth of the pack increases the radius of the ring increases, and interior fills in. Figure 5(a,b) shows the radial distribution of intensity for hcp crystals of $`9`$ and $`25`$ layers in depth in response to a quick impulse force over a small area. Figure 5(c) shows the position of the peaks in intensity as a function of the depth of the pack with the error bars corresponding to the width of the peak at half of the maximum value. For hcp crystals $`25`$ layers in depth the intensity is nearly uniform and thus no lower bound is shown. Additional data points are included corresponding to experimental runs with a quick impulse applied to a large area and to runs with a slowly applied force over a large area. The sloped line indicates an angle of $`27.4`$ degrees which is the average of the angles down the edge and the face of a tetrahedron. The peaks in intensity and force are found to lie just inside of this angle as explained in section IV.

### C Top surface

In addition to measuring the response at the bottom of the crystal packs, we also used the carbon paper technique to measure the response at the top of the pack. Figure 6 shows surface plots for the average intensity of the carbon dots as a function of position at the top surface of the crystal in response to a quick impulse applied to an area two beads in diameter. Both fcc and hcp crystals show a characteristic six-fold spoke pattern, which does not depend significantly on the depth of the crystal. The magnitude of the average intensity in these spokes is approximately ten times less than the intensity in the regions of large force on the bottom surface, but they are still significantly stronger than the surrounding background. The spokes extend all of the way to the edge of the analysis region ($`15`$ bead diameters in radius) without a significant decrease in intensity.

## IV Geometrical Interpretation

### A FCC geometry

Many aspects of the patterns of force or intensity seen in figures 1, 2, 4, and 6 can be explained through force balance in the specific geometry of the crystals. Our fcc crystals were oriented with the $`111`$ plane horizontal, so that they consisted of a stack of triangularly ordered planes. This is equivalent to a typical cannon ball piling. In an fcc crystal the stacking order from plane to plane is such that every third layer lies on top of the first. This results in straight lines of contacts between beads running diagonally down through the pack.

Figure 7 shows the geometry of a piece of a perfect fcc crystal. If a downward force is applied to the bead at the top of this pyramidal pile, then the force is transmitted to the three beads below it in the second layer. Essentially, the bead at the top of the pyramid is supported by a tripod, with three lines of beads stretching downward through the pack. Only the beads in this tripod will feel a force in response to the force on the top bead. All other beads inside and outside of the pyramid are shielded from that force. Thus, at the bottom of the crystal there are only three beads which feel any force. These beads are arranged in a triangular pattern with the size of the triangle proportional to the depth of the pack. The legs of the tripod in the fcc crystal are at an angle of $`35.3`$ degrees with respect to the vertical. Given this angle the position of the beads which feel the applied force is determined. The solid line in figure 3(c) shows the predicted angle of $`35.3`$ degrees.

### B disorder

If there is a small amount of disorder in the crystal, then the positions of the beads within the tripod legs will not be perfectly aligned. Thus, in order to balance the forces vectorally, there must be some force scattered off of the tripod. For a small degree of disorder, the majority of the force will still be transmitted along the tripod legs, with a small amount scattered to one or more additional neighboring beads. Figure 7(c) represents an fcc crystal with one point of disorder in the second layer. Most of the force is transmitted along the tripod legs, but a small amount is scattered along the blue line. If there is no further disorder, then this scattered force will propagate along its own line of beads all of the way to the bottom of the pack, and arrive at the bottom surface along one of the lines connecting the three ends of the tripod legs, as shown in figure 7(d). Using this as a model for the effect of one point of disorder we find that, when considering a small degree of disorder throughout the system, the majority of the force is transmitted along the tripod legs, with a small amount scattered down the sides of the pyramid and arriving at the bottom surface along the edges of the triangle formed by the tripod legs. Figure 7(e) shows the first order effect of scattering within an fcc crystal. The red lines represent strong force, and the blue lines represent weaker forces resulting from first order scattering. The pattern of forces at the bottom of the crystal, shown in 7(f) is seen in the images in figures 1 and 2.

Disorder can also explain the appearance of the regions of large force at angles slightly smaller than $`35.3`$ degrees as seen in figure 3(c) as scattered forces arriving preferentially inside the triangular pattern will skew the regions of large force towards the inside of the pyramid. To first order, no force exists outside of the pyramid.

If at points of disorder the force scatters to beads in plane rather than to the plane below, then without additional scattering the scattered force propagates straight outward to the sides of the container. For nonclose-packed systems, in which beads in-plane are not in contact, it is also possible for the forces to scatter upwards through the pack. This would produce a spoke-like structure, but the length of the spokes would be expected to depend on the depth of the pack which is not seen experimentally.

The second-order effect of disorder on the transmission of forces is depicted in figure 8. Figure 8(a) shows a piece of a three layer fcc crystal. The red lines represent the strong forces along the tripod discussed earlier. The blue lines are the result of first order scattering with both in-plane and downward scattering shown. The green lines represent some of the scattering from the scattered forces, or second order scattering. Of interest is the scattering from the in-plane forces. These forces can scatter upwards, downwards, or in-plane through the pack. Figure 8(b) shows the forces at the top surface of a five layer deep fcc crystal, where second order scattering produces a six fold spoke pattern extending radially outward from the center after a one-bead gap. These spokes extend outward to the edge of the system and do not depend on the depth of the crystal. Similar experimental spoke patterns are shown in figure 6. Figure 8(c) shows the pattern of forces at the bottom surface of the crystal where the effect of second order scattering is to fill in the triangle with very weak forces and extend the edges of the triangle beyond the vertices. For small amounts of disorder the magnitude of the forces from second-order scattering is much smaller than those previously discussed. Nevertheless figure 2(a) shows faint forces along the lines of the triangle edges beyond the vertices.

### C HCP Geometry

The hcp crystals that we studied are also oriented so as to consist of stacked planes of triangular order. Unlike the fcc crystals the stacking of the planes is such that every other plane lies directly above the first. The top two layers of an hcp crystal and fcc crystal are identical, and the force transmission will be the same. However the difference in position of the beads in the subsequent layers causes the force transmission to be markedly different. In an hcp crystal there are not straight lines of beads downward through the pack, so even in a perfect crystal the force cannot be transferred to a single bead and still preserve vector force balance. Furthermore, the nature of the geometry of the hcp crystal is such that the force on a bead in the second layer received from the top bead cannot be transferred to the three beads below. This would require a negative force (see figure 9) or adhesion, which is not present in dry granular systems. The result is that a pack such as that shown in figure 9(a) is not stable to an external force when there is no friction and the weight of the individual grains is neglected. If only the weight of the beads is considered then a perfect hcp packing is stable down to five layers after which it develops the same instability. In order to keep the beads from moving there must be additional support. Inside a larger crystal this instability demands that forces are transmitted to beads in-plane for close-packed crystals and that forces are transmitted upwards for nonclose-packed crystals. If forces are allowed to be transmitted to neighbors in-plane then figure 9(b) shows the propagation of forces for a small piece of a perfect hcp crystal in response to a downwards force at the top. In-plane forces will continue directly to the walls without additional scattering in a perfect crystal. First order scattering will result in a spoke like pattern at the top surface similar to that seen in the fcc crystals. The expected pattern at the bottom surface is a ring of maximum force, whose radius grows with additional layers of depth at angles of $`19.5`$ degrees and $`35.3`$ degrees for alternating layers. The solid black line in figure 5 shows the expected radius of the ring as a function of depth based on an average angle of $`27.4`$ degrees. As the depth of the pack increases the interior of the ring also fills in. This process occurs even in a perfect hcp crystal. Disorder has little effect on the pattern of forces on the bottom surface, but does result in a spoke pattern on the top surface, as is seen experimentally in figure 6.

## V conclusions

We have observed that in the highly ordered regime the response to a localized force is greatly dependent on the structure of the beads within a granular pack. We found distinctly different responses from fcc crystals and hcp crystals, both of which behave differently than amorphous packs. Amorphous packs are found to respond with a large central peak. fcc crystals respond with three regions of large force at the bottom surface, while hcp crystals show rings of maximum force. These structures can be explained by force balance in the geometry of the crystal. These findings are consistent with the two-dimensional results of Geng et al. .

The texture of a granular pack is crucially important for determining the average pattern of stress transmission, as has been suggested by Rajchenbach . The fcc crystal structure, without disorder, leads to a tripod of force chains that is the dominant feature of the force profiles. Small amounts of disorder create secondary structures at the bottom surface and are responsible for a six-fold spoke pattern of forces at the top surface. The hcp crystal structure does not allow for straight line propagation of forces downward through the pack. Even without disorder the forces must split at every layer.

The splitting of forces within an hcp pack can be viewed in analogy to the forcechain splitting model of Bouchaud et al. . In a perfect hcp crystal the force splitting is driven by the orientations of neighboring grains and is thus nearly everywhere uniform, in constrast to the forcechain splitting model in which randomly placed defects drive forces to split into random angles. Nevertheless, the average response of the system to multiple splittings is found to have similarities across these two systems. The forcechain splitting model has been shown to lead to a central minimum of force at the bottom surface in response to a localized force at the top for small packing heights, which is consistent with our observations on hcp crystals of small heights. In two-dimensional simulations of forcechain splitting with a small density of defects the response pattern made a transition towards elastic-like behavior with a broad central peak for depths larger than approximately six mean free paths of the force . We find that for three-dimensional hcp crystals larger than approximately ten layers in depth the central minimum of force fills in leading to a broad central peak response. This experimental result in three dimensions with splitting generated uniformly throughout the pack by contact orientations agrees qualitatively with the predictions of the forcechain splitting model simulations in two dimensions where randomly placed defects generate the splitting.

Finally, it is important to note that previous experiments have shown that the crystal structure does not affect $`P(F)`$, the probability distribution of interparticle contact forces . We conclude that the fluctuations of forces within a single pack are dominated by microscopic details of the granular pack and are not influenced by the arrangement of the beads. Nevertheless, the average response to a localized force is greatly influenced by this structure.

## VI Acknowledgments

We would like to thank Adam Marshall, Daniel Blair, Daniel Mueth, and Milica Medved for their assistance with this project. This work was supported by NSF under Grant No. cts0090490 and by the MRSEC Program of the NSF under Grant No. DMR-9808595.


text:
# Electron transport in a mesoscopic superconducting / ferromagnetic hybrid conductor

## 1 Introduction

The question whether superconductivity can be induced in a ferromagnetic metal is of fundamental and practical importance. At the junction of a Ferromagnetic metal (F) with a Superconductor (S), the superconducting order parameter is predicted to oscillate and decay rapidly in the ferromagnet as the distance to the F/S interface increases. The natural length scale is the exchange length $`L_{exch}=\sqrt{\mathrm{}D/\mu _BH_{exch}}`$, where $`D`$ is the electron diffusion constant and $`H_{exch}`$ the exchange field expressed in Tesla. The latter expression holds in the dirty limit $`L_{exch}>l_e`$, where $`l_e`$ is the elastic diffusion length. Physically, these effects occur because of the wave vector difference between spin-up and spin-down electrons at the Fermi level FFLO ; Beasley . $`L_{exch}`$ is the length over which the two electrons of an Andreev pair get a phase difference of $`\pi `$. Oscillating behaviors were recently detected in measurements of the density of states of a F/S junction Kontos  and in the Josephson supercurrent of a S/F/S junction Ryazanov . These experiments involved F layers with a thickness of the order of the exchange diffusion length $`L_{exch}`$ which was actually made rather large by choosing a ferromagnetic metal with a small exchange field $`H_{exch}`$. In conventional ferromagnetic transition metals (Co which is used here, Ni, Fe, …), the exchange energy $`\mu _BH_{exch}`$ is large and greatly overcomes the thermal energy $`k_BT`$ at cryogenic temperatures. The corresponding exchange diffusion length $`L_{exch}`$ is very small, of the order of a few nanometers.

Many recent experiments involved mesoscopic F/S junctions with a micron-scale ferromagnetic conductor made of Co or Ni. Surprisingly, large proximity effects were observed in the transport properties PetrashovFerro ; Giordano ; Giroud98 ; PrlPetr  of the samples with a transparent interface. When compared to the conventional superconducting proximity effect occurring in non-magnetic metals, JLTPCourtois  this behavior suggests that the relevant length scale is much larger than the expected coherence length $`L_{exch}`$. In the case of interfaces with an intermediate or low transparency, the effect was shown to be restricted to the interface, and was described within an extended BTK model Aumentado . Nevertheless, it was surprising that the fit values of the parameter $`Z`$, which is directly related to the interface transparency, varied so little in comparison with the wide range of interface resistance.

Alternative explanations, ignoring the proximity superconductivity in the ferromagnetic metal, have been proposed. In the case of a transparent interface, the spin accumulation at the F/S junction vanWees ; Falko ; Belzig  could contribute significantly to the anomalous transport properties, together with the Anisotropic Magneto-Resistance (AMR) AMR . Spin accumulation arises near a F/S interface because of the mismatch between the unpolarized pair current in S and spin-polarized single electron current in F. An excess population of minority-spin electrons therefore develops in F in the vicinity of the interface over a length scale set by the spin-flip diffusion length $`L_{sf}`$. Another relevant mechanism could be the competition between the fringe field of the ferromagnet and the diamagnetism of the superconducting electrode. This may result in an inhomogeneous magnetic field distribution and affect locally AMR and/or the Hall effect in the ferromagnet. The AMR stems from spin-orbit coupling in the ferromagnet and results in an anisotropy of the resistivity when the angle between the local magnetization and the current flow changes. This enables the observation of magnetization reversal processes in the resistance of small magnetic particles APLAumentado .

The microscopic mechanism inducing superconductivity in F close to the S interface is the Andreev reflection where an incident electron is reflected into a phase-correlated hole of the same spin. This is equivalent to creating an Andreev pair of electrons with opposite spins. Obviously, this process will be affected by the spin polarization of the ferromagnetic metal and will even disappear in the case of a fully-polarized metal Beenakker ; Soulen ; Upadhyaya . In this respect, an open question is the role of the magnetic domain structure. Crossed Andreev reflections Feinberg  may appear at the F-S interface close to a domain wall separating in F two magnetic domains with opposite magnetization Melin . In the case of a ferromagnet with an inhomogeneous magnetization, it was also proposed that the spin-triplet component of the superconducting wave function can have a strong amplitude Volkov ; Kadigrobov . Interestingly, this component should exhibit a slower spatial decay than the usual singlet component. This inhomogeneous magnetization hypothesis is relevant because of both the shape anisotropy of micro-fabricated structures and the effect of the diamagnetism of the S electrode.

These open questions show the need for further investigation of transport in a mesoscopic ferromagnetic conductor connected to a superconductor. Here, we report on transport measurements of submicron Co-Al hybrid structures. Compared to our previous experiments Giroud98 , we modified both sample dimensions and geometry to focus on transport properties in the ferromagnet itself, near the superconducting contact. We observed a large resistance drop in samples with a highly transparent interface. We studied the dependence of this effect with the temperature, the magnetic field and the bias current. Our main result is the observation of a resistance asymmetry between the two nominally identical branches of the Co wire.

## 2 Samples description

### 2.1 Fabrication

The samples geometry (see Fig. 1) was designed for measuring the ferromagnet resistance in proximity to the superconducting contact, with zero net current through the interface. The Al contact is deposited on a lateral Co ”finger”, rather than directly on top of the Co wire itself, in order to minimize spurious current density redistribution effects when Al becomes superconducting NSa2D ; Belzig . As there are evidences that the superconducting contact influences transport only in the vicinity of the interface, the Co strip length was chosen as short as possible, namely 400 nm.

The sample fabrication process was chosen as to keep F/S contact resistances as small as possible. We used a two-step lift-off process with in-situ Ar ion etch. Co was deposited first on the silicon substrate, in order to avoid step edges which could modify magnetization anisotropy and pin magnetic domain walls. A 50 nm layer of Co was e-beam evaporated at room temperature through a PMMA mask in a base vacuum below $`10^7`$ mbar. An in-situ Ar ion milling of the Co surface was performed just before the 100 nm Al layer evaporation through the second PMMA mask.

### 2.2 Characterization

The deposition conditions together with the Co thickness are expected to result in an in-plane magnetization, which was confirmed by the AMR data. This orientation makes it easier to induce magnetization reversal under an applied magnetic field. The expected typical domain size is of the order of 100 nm, i.e. roughly comparable to the wire width. The Co resistivity was reproducibly high, in the 80 $`\mu \mathrm{\Omega }`$.cm range, whereas Al residual resistivity did not exceed 2 $`\mu \mathrm{\Omega }`$.cm. This corresponds to an electron diffusive mean free path $`l_e`$ of 1 nm in Co and 20 nm in Al. This also gives an estimated coherence length for superconducting correlations $`L_{exch}`$ in Co of 3 nm. The Al superconducting coherence length and London penetration depth are of the order of 0.12 $`\mu `$m and 0.18 $`\mu `$m respectively. Overall, the Al resistive transition at $`T_c`$ 1.3 K did not seem strongly affected by the proximity of the Co wire, except for a depressed critical current, roughly 3 times lower than for pure Al strips or microbridges of comparable dimensions.

We are mainly interested by the case of low resistance interfaces, as it is presumably the only case where the superconductor has a strong influence on the ferromagnet. Interface resistance can be probed using two probes on the Al electrode and two probes at both ends of the Co strip. Due to the sample geometry, we can only measure a global resistance $`R_j`$ which is the interface resistance in series with the Co lateral finger resistance and the Co spreading resistance. In the following, we consider experimental data collected on two samples (A and B) with a low resistance interface. The distance d between the main Co wire and the Al contact, i.e. the length of the lateral Co finger not covered by Al, was estimated from SEM micrographs to d $``$ 50 nm in sample A and d $``$ 100 nm in sample B. The global ”junction” resistance $`R_j`$ was about 17 $`\mathrm{\Omega }`$ in sample B. From the sample dimensions (see Fig. 1) and Co resistivity, we estimate that the Co finger itself is the main contribution. We can thus deduce that our interface specific resistance is below 6 $`m\mathrm{\Omega }.\mu `$m2 in these samples with a transparent interface. For comparison, we also discuss data from one sample (C) with a degraded interface. The resistance $`R_j`$ is 103 $`\mathrm{\Omega }`$, which leads to an estimated interface specific resistance of 600 $`m\mathrm{\Omega }.\mu `$m2.

## 3 Experimental results

### 3.1 Measurement procedure

We studied the electron transport properties of several samples down to a temperature of 30 mK. A magnetic field was applied in the sample plane, either parallel or perpendicular to the current in the Co wire. An experimental run with a field perpendicular to the substrate plane was also performed on some samples, to check that the Co AMR had the behavior expected for in-plane magnetization. In all our measurements, a d.c. current bias plus an a.c. current modulation below 100 nA were applied through the same contacts. Low-pass filters were inserted on the cryostat feedthroughs, as well as a d.c. rejection filter at the input of the lock-in amplifier. Depending on the wiring (see Fig. 1 right part), we investigated the resistance $`R_{Co}`$ of the Co wire between the two reservoirs, that of one of the halves of this wire $`R_+`$ or $`R_{}`$, or the junction resistance $`R_j`$.

### 3.2 Temperature dependence

Fig. 2 shows the temperature dependence of the Co conductor in samples A, B and C. As the temperature is decreased below the critical temperature of Al, the resistance of sample A first drops by about 1 $`\%`$. It drops again by about 10 $`\%`$ below 0.2 K. Qualitatively similar results were obtained in other samples with a low resistance Co-Al interface. This is the case of sample B, which was patterned on the same wafer than sample A, but then the temperature of the large resistance drop is much higher whereas the first resistance drop just below 1.3 K is not so clearly visible. Therefore the characteristic temperature of the resistance drop, as well as the shape of the curve, appear to be sample-dependent. We do not yet know precisely which factors monitor this variation. Let us point out that the total resistance drop magnitude is nearly 12 $`\%`$ in samples A and B. This effect is not observed at large bias current.

In the case of sample C with a poor interface transparency, the superconducting transition has almost no effect on the ferromagnet resistance $`R_{Co}`$, as shown in Fig. 2. It only weakly affects the junction resistance $`R_j`$, which first decreases of about 0.4 $`\%`$, and then increases on cooling down (see Fig. 3).

### 3.3 Magnetoresistance

Fig. 4 shows the magneto-resistance of sample A at low bias current for a magnetic field applied in-plane along the Co wire, i.e. parallel to the current path. At low field ($`H<`$ 150 mT), the magnetoresistance shows a small amplitude (less than 1 $`\%`$) and a significant hysteresis. These two features suggest that this low-field behavior is due to Anisotropic Magneto-Resistance (AMR). The demagnetized multi-domain state reproducibly exhibits a larger resistance than the higher field state. The relatively sharp jumps at $`\pm `$ 70 mT are the signature of the cobalt coercitive field $`H_{coer.}`$. In the high-field regime, we observe a large resistance increase up to the normal state residual resistance. This increase is the signature of the Al transition to the normal state at the Al superconducting critical field $`H_{cS}`$ 210 mT.

Measurements on sample C (not shown) showed a small 0.6 $`\%`$ positive magneto-resistance when field and current are in-plane but perpendicular to each other. This anisotropy, as well as the amplitude of the resistance jumps, are consistent with the effect of the Co AMR in the case of an in-plane magnetization. The small resistance increase at very low temperature for sample C (see Fig. 3) can therefore be explained by a modification of the Co AMR induced by Al diamagnetic shielding.

Let us note a significant difference with our previous experiments Giroud98  on 2 $`\mu `$m-long samples. In these long samples, magnetoresistance showed a smooth variation over a relatively broad field range. In the present short ”T-shaped” samples, we reproducibly observe very well defined jumps at $`\pm H_{coer.}`$. This strongly suggests that our short wires contain a small number of domain walls which depin at this field, whereas long samples contain a larger number of domains walls, which will rotate or move along the sample in a broader field range.

### 3.4 Differential resistance

The differential resistance was investigated by superposing a d.c. current bias current to the a.c. modulation. Fig. 5 displays the total differential resistance measured in a standard four probes configuration between the ”+” and the ”-” ends of the Co wire. In sample A, we observe a resistance dip for current bias below 0.2 $`\mu `$A which mimics the resistance drop observed below 0.2 K. At this point, the voltage across sample A is about 20 $`\mu `$V, roughly a factor 10 below the expected Al gap. The high-bias resistance value coincides with the value obtained at high temperature, above the critical temperature of Al. The relative amplitude of the variation is in the 10 $`\%`$ range. A similar behavior was reproducibly observed in sample B and several other samples. The current range and the profile of the resistance dip appear to scale with the characteristic temperature of the resistance drop. No variation of the junction differential resistance was observed at low bias except for temperatures close to $`T_c`$.

Because of the relatively high resistivity of our Co film, we have to be careful about heating effects. At the lowest temperature (30 mK), we expect that the heat flow is essentially evacuated along the Co strip. At a temperature of 0.1 K, the thermal conductance of the Co strip, Co/Si interface, and Al strip, should be respectively about 100, 50 and 10 pW/K. The Joule power dissipated by the Co strip is only 4 pW at the 0.2 $`\mu `$A current bias required to suppress the resistance decrease. This means that at this bias the sample A cannot be heated above 0.1 K. We conclude that Joule heating is not sufficient to explain the differential resistance variation.

We also applied a $`\pm `$ 124 mT field along the Co strip, parallel to the current path. The value $`\pm `$ 124 mT was chosen as to be above the Co coercitive field $`H_{coer.}`$ 70 mT and below the Al superconducting critical field $`H_{cS}`$ 210 mT. The Co conductor was therefore presumably close to magnetization saturation, although it may not yet be single domain. The current dependence of $`R_{Co}`$ shown in Fig. 5 only shows a small difference compared to the zero field case. This difference is of the same order than the Co magnetoresistance jump at $`H_{coer.}`$, i.e. close to 1$`\%`$.

### 3.5 Resistance asymmetry

We now come to our main result. Figure 6 top part shows sample A differential resistances of the left and right parts of the Co wire ($`R_+`$ and $`R_{}`$), as well as the total resistance $`R_{Co}=R_{}+R_+`$ and the resistance difference $`R_{}R_+`$. The current path between the two ends of the Co strip remained unchanged, so that the net current through the interface is always zero. $`R_+`$ and $`R_{}`$ represents the differential resistance measured respectively with voltage probes at the ”+” Co end and the Al contact, or between the Al contact and the ”-” Co end (see Fig. 1).

We observe a non-trivial asymmetry between $`R_+`$ and $`R_{}`$ in the low bias regime. In the current range where Co resistance is depressed by the superconducting contact, $`R_{}`$ is lower and $`R_+`$ higher than its respective high-bias value. The relative variations $`\mathrm{\Delta }R_{}/R_{}`$ and $`\mathrm{\Delta }R_+/R_+`$ are significantly stronger (up to 35 $`\%`$) than the total Co wire variation $`\mathrm{\Delta }R_{Co}/R_{Co}`$ (about 12 $`\%`$). At high bias, the difference between the values of $`R_+`$ and $`R_{}`$ may be explained by the inhomogeneities in the sample cristalline micro-structure. Let us note that the differential resistance remains symmetric with respect to the d.c. current bias, i.e. a current reversal produces a voltage sign reversal. Within a small experimental error, the sum of $`R_+`$ and $`R_{}`$ is always equal to the total differential resistance $`R_{Co}`$. Similar results were reproducibly observed on several samples, including sample B (Fig. 6 bottom).

## 4 Discussion of the results

First, let us comment on a difference between these results and our previous experiments Giroud98 : we no longer observed a re-entrance of the metallic resistance. Here, the estimated Thouless energy of the Co conductor is $`ϵ_c=\mathrm{}D/L^2=0.1K`$, and a resistance minimum should have been observed close to 0.5 K, which is not the case. This contradicts the interpretation we proposed earlier, namely the occurrence of a re-entrant proximity effect in the resistance, similarly to the non-ferromagnetic metal case JLTPCourtois . We come to the conclusion that the resistance minimum in our previous experiments may rather result from the competition between two opposite mechanisms: a resistance drop induced by the superconducting contact which we observe much more clearly in our shorter new samples, and a resistance upturn of different origin. Spin accumulation effects constitute an obvious candidate Belzig , but Anisotropic Magneto Resistance (AMR) cannot be excluded.

The first question we have to care about is whether the resistance drop is indeed occurring in the ferromagnet, or is merely a consequence of current redistribution in the superconducting short-circuit NSa2D . It is important to note that in our samples, the sheet resistance of Al above $`T_c`$ is only 0.2 $`\mathrm{\Omega }`$ per square, which is much smaller than the sheet resistance of Co ($``$ 20 $`\mathrm{\Omega }`$ per square). Even normal Al acts as a shunt. One can model our two-dimensionnal sample with an array of resistances of the order of the Co and Al sheet resistances. The result of this analysis is that the resistance drop expected at Al superconducting transition should not exceed the Al sheet resistance (0.2 $`\mathrm{\Omega }`$). This is clearly much smaller than the experimental resistance drop, which exceeds 10 $`\mathrm{\Omega }`$. It is also known that the AMR is of the order of only 1 to 2$`\%`$ in ferromagnetic 3d transition metals such as cobalt. We indeed observe this AMR effect in the low-field magnetoresistance (Fig. 4). Thus we are confident that the resistance drop observed in Co, about 10 $`\%`$ above the residual normal state case, is neither a trivial current redistribution effect in the superconducting short-circuit nor a simple AMR effect. The behavior of otherwise identical samples with respectively low or high junction resistance is so clearly distinct that we can conclude that interface transparency plays a key role. The comparison between Fig. 2 and 3 even shows that the sign of the resistance variation can be reversed, as was observed in Ref. PrlPetr . In the case of a weakly transparent interface, our results are also compatible with the conclusion of Ref. Aumentado  that no proximity effect appears in the bulk of the ferromagnet. A significant influence from the superconductor on the ferromagnet occurs when the interface is transparent enough, and only in this case.

The differential resistance asymmetry when changing the voltage probes configuration is somewhat surprising since the sample were fabricated as symmetric. We have to consider physical phenomena which may contribute differently to the resistances $`R_+`$ and $`R_{}`$ of the two sample halves. From the sample geometry shown in Fig. 1, we see that $`R_+`$ and $`R_{}`$ may include a Hall effect contribution, related to the local field in the ferromagnet or to its magnetization (anomalous Hall effect) over the Co transverse dimensions (wire width plus finger length). The Hall voltage will contribute with opposite sign to $`R_+`$ and $`R_{}`$, but should not contribute to the total resistance $`R_{Co}`$, as in this case voltage probes are aligned along the current lines. A d.c. current bias as low as 0.2 $`\mu `$A cannot significantly affect the magnetization and the domains structure tsoi ; wegrowe01 , so that the Hall resistance should remain constant in the experimental bias range. Moreover, we observed that samples with a high resistance interface do not exhibit such a significant resistance asymmetry, although their magnetization and coercitive field are not significantly different. Therefore, the differential resistances and their asymmetry cannot be explained by the classical anomalous Hall effect in the ferromagnet, but are related to superconductivity.

In the interface region, electrons diffuse between Co and Al. In the case of a positive d.c. bias electrons coming from the ”+” side have a higher energy than those going on the ”-” side. An out-of-equilibrium energy distribution will develop near the interface even with a zero net current JLTPCourtois . Because of the mismatch between spin-polarized current in F and unpolarized pair current in S, any current between S and F will be also associated to a non-equilibrium spin polarization, in a distance range from the interface determined by the spin flip diffusion length. This spin polarization is energy-dependent, should be maximum at zero bias and vanish above the superconducting gap. With this simple picture in mind, one might understand that $`R_+`$(I) and $`R_{}`$(I) may be different at low bias. Nevertheless, we expect that $`R_{}`$(-I) should behave as $`R_+`$(I) if the sample is symmetric. This is clearly not the case, one exhibiting a maximum at zero bias and the other one a minimum.

Therefore, the difference between $`R_+`$ and $`R_{}`$ must stem from a physical asymmetry in the sample itself, the most obviously possible one being the magnetic domains structure in an otherwise symmetric geometry. Since the T shape results in a complicated shape anisotropy, it is very likely that a non-symmetric magnetic domain structure is present in the central region of our sample. The amplitude and the sign of this asymmetry will be of course sample-dependent. If we assume that there is only a small number of domains in our short samples, and if the spin flip diffusion length is not much smaller than domain size (e.g. in the 200 nm range), the chemical potential drop may be different for electrons travelling from one Co contact (or the other) to the Al electrode, depending on the magnetic domain structure in the $`R_+`$ or $`R_{}`$ section of the cobalt wire. Nevertheless, the actual effect of the superconductivity on the electron transport in Co and the physical origin of the resistance drop remains undetermined. It could be a long-range superconducting proximity effect like the predicted long-range triplet component Volkov ; Kadigrobov , spin accumulation effect in close relation with the sample geometry Falko ; Belzig , or Andreev reflections of electrons of opposite spins in adjacent ferromagnetic domains of different magnetization Feinberg ; Melin .

## 5 Conclusion

In this work, we brought new experimental evidence for large resistance decrease in hybrid Ferromagnetic / Superconducting devices. This effect is clearly distinct from the Anisotropic MagnetoResistance (AMR) of smaller amplitude. A high interface transparency was found to be necessary for observing large effects. We suggest the relevance of the magnetic domain structure in the transport properties in the vicinity of a F/S contact. Further work on transport properties of F/S junctions with a high control of the magnetic domain structure is in progress .

We thank W. Belzig, R. Mélin and M. Viret for discussions and C. Veauvy for help in cryogeny. We acknowledge discussions within the ”Dynamics of Superconducting Nanocircuits” TMR network and the ”Mesoscopic Electronics” COST P5 action. This work is supported by the French Ministry of Education and Research under an ACI contract.