text:
# Electronic and structural properties of superconducting MgB2, CaSi2 and related compounds

## I Introduction

The recent discovery of superconductivity at $`39`$K in the simple intermetallic compound MgB2 is particularly interesting for many reasons. In the first instance, this critical temperature is by far the highest if we exclude oxides and C60 based materials. Furthermore, B isotope effect suggests that MgB2 is a BCS phonon-mediated superconductor, with $`T_c`$ above the commonly accepted limits for phonon-assisted superconductivity. Another reason of interest is that the AlB2-type structure of MgB2 (which can be viewed as an intercalated graphite structure with full occupation of interstitial sites centered in hexagonal prisms made of B atoms) is shared by a large class of compounds (more than one hundred), where the B site can be occupied by Si or Ge, but also by Ni, Ga, Cu, Ag, Au, Zn and Cd, and Mg can be substituted by divalent or trivalent $`sp`$ metals, transition elements or rare earths. Such a variety is therefore particularly suggestive of a systematic study of diborides and structurarly related compounds. To confirm this interest, a recent study showed that superconductivity is supressed by addition of more than about 30% of Al atoms. Pressure effects also lead to a depression of superconductivity.

Early studies of the electronic structure of diborides by the OPW and tight-binding methods have shown similarity with graphite bands. A more systematic study of $`\mathrm{AlB}_2`$–type compounds by the full-potential LAPW method showed in simple metal diborides the presence of interlayer states, similar to those previously found in graphite, in addition to graphite–like $`\sigma `$ and $`\pi `$ bands. Transition metal diborides, on the other hand, show a more complex band structure due to the metal $`d`$ states. A very recent electronic structure calculation on MgB2 shows that metallicity is due to B states. Phonon frequencies and electron-phonon couplings (at the $`\mathrm{\Gamma }`$ point) lead these authors to suggest phonon-based superconductivity. Further, new calculations on related compounds ($`\mathrm{LiB}_2`$ and graphitic-like B) addressed the peculiar contributions to conductivity coming from the interplay of 2D and 3D Fermi surface features. An alternative mechanism has been proposed by Hirsch, based on dressed hole carriers in the B planes, in the presence of almost filled bands.

Recently, superconductivity at $`T_c14`$K was observed in $`\mathrm{CaSi}_2`$. In this compound, high pressure has been effectively used to tune both structural and physical properties. In particular, Bordet and co-workers  identified a new structural phase transition for $`\mathrm{CaSi}_2`$ from trigonal to hexagonal $`\mathrm{AlB}_2`$-type structure, with superconductivity at $`T_c=14K`$, one of the highest transition temperatures known for a silicon based compound. The occurrence of this structural phase transition indicates that Si shifts from an $`sp^3`$ to an $`sp^2`$-like hybridization, and experiment points to a clear enhancement of $`T_c`$ with high pressure. Many electronic structure studies have been carried out for the different phases of disilicides. $`\mathrm{CaSi}_2`$ in the $`\mathrm{AlB}_2`$ structure was studied by Kusakabe et al., but a final answer on the stability of this phase under pressure is still lacking.

As previously pointed out, the $`\mathrm{AlB}_2`$ family of compounds provides a wide playground for the search of new superconductors, and we intend to provide here a contribution, based on a detailed study of electronic and structural properties together with phonon spectra, for several diborides and for $`\mathrm{CaSi}_2`$, in the trigonal and hexagonal structures. The questions that we want to address are the following: (i) how do the electronic properties of MX2 depend on the atomic species M; namely, can rigid band schemes provide the correct trend, or is the metal effect mostly mediated by lattice parameter changes (chemical pressure). (ii) Ref. shows the presence of both B $`\pi `$ Fermi surfaces (FS) giving rise to 3-dimensional states, and of small 2-dimensional B $`\sigma `$-bonding pockets. It is of course interesting to know which of the FS is mostly involved in superconductivity. The trends observed in the band structures, together with experimental studies, will contribute to elucidate this point. In the same way, changes in the Fermi level density of states with chemical composition will possibly correlate with $`T_c`$ values. (iii) Pure CaSi2 shows a trigonal instability, joined by a buckling of Si2 planes, which is apparently strongly reduced under high pressure, or by Be substituting some of the Si sites. The presence, or the proximity, of structural instabilities (investigated through phonon spectra), will hopefully provide hints to further experimental search for new materials.

In this paper, we present calculations for MgB2, $`\mathrm{Mg}_{0.5}\mathrm{Al}_{0.5}\mathrm{B}_2`$, BeB2, CaSi2 and $`\mathrm{CaBe}_\mathrm{x}\mathrm{Si}_{2\mathrm{x}}`$ ($`x=1`$), based on the local density approximation (LDA) to the density functional theory. We use both the Full Linearized Augmented Plane Wave method, and the norm-conserving pseudopotential method. In Sect. II we give some computational detail; in Sect. III and IV we discuss the electronic properties for all the compounds considered; we finally draw our conclusion in Sect. V.

## II Computational and structural details

Our calculations were performed using the “all–electron” full potential linearized augmented plane waves (FLAPW)  method, in the local density approximation (LDA) to the density functional theory. In the interstitial region we used plane waves with wave vector up to $`K_{max}`$ = 3.9 a.u. and 4.1 a.u. for disilicides and diborides, respectively (we used the larger $`K_{max}`$ = 4.4 a.u. for BeB2 because of the smaller spheres imposed by the smaller lattice parameters). Inside the muffin tin spheres, we used an angular momentum expansion up to $`l_{max}`$ = 6 for the potential and charge density and $`l_{max}`$ = 8 for the wave functions. The Brillouin zone sampling was performed using the special $`k`$–points technique according to the Monkhorst-Pack scheme, and also the linear tetrahedron method with up to 120 $`k`$–points in the irreducible Brillouin zone. We used muffin tin radii $`R_{MT}`$ = 2.2, 1.8, 1.59 and 2.1 a.u., for Ca, Si, B, and the remaining elements respectively (we used $`R_{Be}`$ = 2 a.u. and $`R_B`$ = 1.5 a.u. in BeB2).

The crystallographic structure of $`\mathrm{AlB}_2`$ (C32) is hexagonal, space group $`P6mmm`$ with graphitic B planes, and 12-fold B–coordinated Al atoms sitting at the center of B hexagonal prisms. The unit cell contains one formula unit, with Al atom at the origin of the coordinates and two B atoms in the positions ($`𝐝_1=(a/2,a/2\sqrt{3},c/2)`$, $`𝐝_2=𝐝_1`$). The trigonal distortion in CaSi2 corresponds to a buckling of Si planes leading to a $`P3\overline{m}1`$ space group symmetry), with one internal structural parameter $`z`$ (the internal coordinate of B atoms).

To calculate plasma frequencies, Fermi velocities, Hall coefficients, and the densities of states (DOS) near $`E_F`$, we have performed a spline fit of the ab-initio energy bands. The fitted bands have then been computed over 1283 points in the reciprocal space unit cell. The linear tetrahedron method has been used to compute DOS and the other quantities.

Phonon frequencies were computed using the ABINIT code  by means of Linear Response Theory  (LRT) approach in the framework of pseudopotentials plane wave method. Troullier and Martins soft norm-conserving pseudopotentials were used for all the elements considered, Mg and Ca $`p`$ semicore states were not included in the valence, their contribution to the exchange and correlation potential was accounted by the non-linear core corrections . The local density approximation to the exchange-correlation functional is used. Plane wave cutoffs of 44 Ry and 30 Ry were found sufficient to converge structural and electronic properties of the diborides and silicides respectively. Brillouin zone integration was performed using a $`20\times 20\times 20`$ Monkhorst-Pack mesh for the structural optimization and a $`12\times 12\times 12`$ mesh was used in the LRT selfconsistent calculation of phonon frequencies. To improve $`k`$–point sampling convergence a gaussian smearing with an electronic temperature of 0.005 Ha was used. The structure used in the LRT calculations was obtained optimizing the lattice structure up to a residual stress of 1 MPa.

## III Diborides

### A MgB2 and BeB2.

We shall start our discussion of electronic properties of diborides by presenting in Fig. 1 the energy bands of the superconducting compound MgB2, which will be the reference for all our further studies (the $`M\mathrm{\Gamma }K`$ lines are in the basal plane, while the $`LAH`$ lines are the corresponding ones on the top plane at $`k_z=\pi /c`$). We used the experimental lattice parameters $`a=3.083`$ Å, $`c=3.52`$ Å. Our results, which are in excellent agreement with those of Kortus et. al., show strong similarities with the band structures of most simple-metal diborides, computed by the FLAPW method. They show strongly bonded $`sp^2`$ hybrids laying in the horizontal hexagonal planes, forming the three lowest ($`\sigma `$–bonding) bands, as the main structure in the valence region. The corresponding antibonding combinations are located around 6 eV above the Fermi level ($`E_F`$).

The bonding along the vertical direction is provided by the $`\pi `$ bands (at $``$ -3 eV and $``$ 2 at $`\mathrm{\Gamma }`$), forming the double bell-shaped structures. Similarly to $`\mathrm{AlB}_2`$-type diborides, and unlike graphite, we notice a large dispersion of the $`\pi `$ bands along the $`k_z`$ direction ($`\mathrm{\Gamma }A`$). This can be understood on a tight-binding approach: the phases of the Bloch functions lead the B $`p_z`$ orbitals on adjacent layers to be antibonding at $`\mathrm{\Gamma }`$, and bonding at $`A`$. The downwards dispersion of the $`\pi `$-bonding band (from $`3`$ to $`7`$ eV along $`\mathrm{\Gamma }A`$) is joined by an opposite upwards dispersion of an empty band, (at $`2`$ eV at the $`\mathrm{\Gamma }`$ point). The analysis of the wavefunctions for this empty band reveals that this band has an interstitial character, and is very similar to those found in graphite, in diborides, in CaGa2, but also in CaSi2. A general description of the bonding can be given by the partial density of states (PDOS), shown in Fig. 2. The B PDOS, in particular, shows bonding and antibonding structures for the $`s`$ and $`p`$ states. The $`p_z`$ states give rise to a rather wide structure responsible for less than half of the DOS at $`E_F`$, which therefore results having predominantly B $`p\sigma `$ character.

In the search of similar materials having the desired superconducting properties, BeB2 is the first natural candidate, as the band filling is expected to be pretty similar, and lighter Be atoms may help providing larger phonon frequencies while keeping similar electronic properties. The first experimental report on BeB2 gives the average lattice constants. These experimental parameters have been used for all the previous calculations of the BeB2. We have optimized the lattice constants, imposing the $`\mathrm{AlB}_2`$ structure, obtaining $`a=2.87`$Å, $`c=2.85`$Å, smaller than in MgB2 and similar to the values estimated by the average experimental values. In particular, we notice that while $`a`$ reduces only by about 7.4%, $`c`$ reduces by about 24% when Mg is substituted by Be, a result which can be understood pretty easily. On one hand, in fact, the optimization of $`\sigma `$-bonds prevents a too drastic variation of the $`a`$ parameter; on the other hand, $`c`$ can change more easily as the vertical bonding, provided by the $`\pi `$ bands, has a large contribution from the metal orbitals. The band structure of BeB2, shown in Fig 3, are pretty similar to those computed by the FLAPW method at the experimentally estimated lattice parameters. We can remark the strong similarities with the MgB2 compound, which includes the presence of $`\sigma `$-bonding hole pockets along $`\mathrm{\Gamma }A`$, but also a few relevant differences. First of all, the shorter lattice parameters lead to wider valence bands, and also to more dispersed $`\sigma `$ bands, in particular the cylindrical hole pocket along $`\mathrm{\Gamma }A`$, which might be relevant for superconductivity. Also, the different energy location of the metal $`s`$–free electron band (lower in MgB2), and the different shape of the $`\pi `$-bonding band (related by the different $`c`$ values) causes a different occupation of this band, especially along $`\mathrm{\Gamma }M`$.

To further investigate these two materials, and their relation to isostructural diborides, we studied the zone center phonon frequencies of MgB2, BeB2, and $`\mathrm{AlB}_2`$, using the linear response approach, with the pseudopotential method. The calculations have been carried out at the theoretical lattice constants (also listed in Table I), as obtained using the same method. The two lowest independent frequencies correspond to a motion of the intercalated metal atoms relative to the rigid B networks: one vertical mode ($`A_{2u}`$), and one doubly degenerate in-plane mode, ($`E_{1u}`$). The two higher frequencies represent the internal motion of the B network itself: again one vertical mode ($`B_{1g}`$), now representing a trigonal distortion producing a buckling of the planes similar to that found in CaSi2, and one degenerate ($`E_{2g}`$) mode giving an in-plane stretching of B-B bonds. Our calculations for MgB2 are in good with other results . The comparison between the MgB2 and BeB2 shows a much higher $`E_{2g}`$ frequency in BeB2. This result cannot be explained through the different atomic masses, as it involves only B motion. It is therefore related to the different coupling of the electrons with the nuclear system. Surprisingly, two of the remaining frequencies are lower in BeB2, namely the $`A_{2u}`$ mode which involves the motion of lighter Be atoms against the B network and the $`B_{1g}`$ mode. This indicates a strong dependence of the phonon spectra of these materials on the electronic structure details. The lowest $`E_{1u}`$ mode has comparable frequencies. If we now go to $`\mathrm{AlB}_2`$, we see again that the $`E_{2g}`$ mode is much higher than in MgB2. Therefore, this mode has his minimum value in the superconducting material. The modes involving the relative motion of the B2 planes as a whole have frequencies comparable to the other materials; the $`B_{1g}`$, on the other hand, is smaller than in MgB2 and BeB2, which is probably a consequence of the larger band filling (see the discussion on CaSi2 below).

We have computed, within a rigid band scheme and using the scheme described, e.g., in Ref., the Fermi velocities, plasma frequencies (along the principal axes of the crystal), and the independent components of the Hall tensor for BeB2 and MgB2. Because of hexagonal symmetry, the non-zero components will be $`R_{xyz}`$, corresponding to the magnetic field along the $`z`$ axis and transport in plane, and $`R_{zxy}=R_{yzx}`$ corresponding to in–plane magnetic field. The results, plotted in Fig. 4 as a function of the hole or electron-type doping, show relatively smooth variations over a wide range of doping, and are nearly coinciding, at zero doping, with the corresponding calculations by Kortus et al.. After comparing the two materials, the following remarks can be done: (i) the density of states at $`E_F`$ is smaller in BeB2. (ii) MgB2 is nearly isotropic in terms of plasma frequencies and Fermi velocities, while BeB2 is expected to show important anisotropy in resistivity, and, in particular, the in-plane conductivity is expected to be smaller. A band-by-band decomposition explains these results in term of very similar (and very anisotropic) contributions in the two compounds from $`\sigma `$ states, and quite different contributions from the $`\pi `$ bands, nearly isotropic in MgB2 and favoring conductivity along the $`c`$-axis in BeB2. These results point out that despite structural similarity, the electronic properties of $`\mathrm{AlB}_2`$-type structures differ substantially from graphite and intercalation graphite compounds. (iii) The Hall coefficients are positive (hole-like) when the magnetic field is along the $`z`$ axis (transport in–plane), while is negative for MgB2 and almost zero but negative for BeB2, when the magnetic field is in–plane. A detailed analysis shows again comparable contributions in the two compounds from $`\sigma `$ bands, while the sign of $`R_{zxy}`$ and $`R_{yzx}`$ comes from a balance of the $`\pi `$ bands. The Hall coefficient of MgB2 has been recently measured by Kang et al.. $`R_H`$, which will correspond in this case to an average of the tensor components, turns out to be positive and to decrease with temperature. At $`T=100`$K they give $`R_H=4.1\times 10^{11}\mathrm{m}^3/\mathrm{C}`$, and an extrapolation of their results at $`T0`$, gives $`R_H6.5\times 10^{11}\mathrm{m}^3/\mathrm{C}`$. If we average, at zero doping, our tensor components, the positive contributions overcome the negative ones, resulting in an average value $`R_H2\times 10^{11}\mathrm{m}^3/\mathrm{C}`$. This value is thus of the correct order, but single crystal measurements are probably necessary to have a significant comparison, avoiding an average of quantities with different sign. It is worth noticing that, if the sign of carriers has a role on superconductivity, BeB2 also has a positive Hall coefficient.

### B Al doping in MgB2.

A recent experimental report by Slusky et al. shows that Al doping destroys bulk superconductivity when the Al content $`x`$ is greater than $`0.3`$. Stimulated by these results, we studied the chemical substitution of Mg with Al. The substitution corresponding to a 50% concentration of Al ($`x=0.5`$) leads to a reduction of the cell parameters, from $`a=3.083`$ Å, $`c=3.521`$ Å for MgB2, to $`a=3.047`$ Å, $`c=3.366`$ Å for $`\mathrm{Mg}_{0.5}\mathrm{Al}_{0.5}\mathrm{B}_2`$. These differences again indicate hardly compressible B-B $`\sigma `$ bonds, and interlayer distances quite sensitive to the intercalated cation. The study of Slusky et al. shows consistent results, and indicates, furthermore, that after a two-phases region for Al substitutions from 10 to 25%, $`c`$ collapses and bulk superconductivity disappears. We first studied this problem in a rigid band approach, bringing MgB2 to the same lattice parameters of $`\mathrm{Al}_{0.3}\mathrm{Mg}_{0.7}\mathrm{B}_2`$. The only apparent change at that electron concentration is the complete filling of the $`\pi `$ bonding band at the $`M`$ point. The $`\sigma `$ bonding bands need a larger electron addition to be completely filled in a rigid band model.

To investigate further this problem, we studied the system with 50% Al concentration. We used in our calculations the experimental values, and in a first approach simulated the disorder (certainly present in the real compound) using an orthorhombic supercell. The corresponding unit cell contains two formula units, and is described by the two orthogonal in plane basis vectors $`𝐚_1=(a,0,0)`$ and $`𝐚_2=(0,a\sqrt{3},0)`$. We started by calculating the energy bands using the MgB2 lattice constants, and then we used the experimental $`\mathrm{Mg}_{0.5}\mathrm{Al}_{0.5}\mathrm{B}_2`$ constants. The results of these calculations are shown in Fig. 5 together with the bands of MgB2, folded into the Brillouin zone (BZ) of the orthorhombic structure. The folding is such that the $`\mathrm{\Gamma }`$ point of the orthorhombic BZ corresponds to the $`\mathrm{\Gamma }`$ and $`M`$ points of the hexagonal BZ, while $`Z`$ corresponds to both the hexagonal $`A`$ and $`L`$. Therefore, the $`\mathrm{\Gamma }Z`$ line collects the $`\mathrm{\Gamma }A`$ and the $`ML`$ hexagonal lines, and we notice immediately the three $`\pi `$ bands (two bonding and one anti-bonding, coming from folding). The $`\mathrm{\Gamma }X`$ line contains the (folded) hexagonal $`\mathrm{\Gamma }M`$ line. If we first look at the frozen structure calculation Fig. 5(c), we see quite similar band dispersions with the obvious splittings, related to the symmetry lowering induced by the different cations, and a different band filling. We notice, however, a lowering of the folded $`\pi `$ antibonding band, related to the different cations: the $`\pi `$ band is in fact pushed downwards by the interaction with the free-electron metal-$`s`$ band, which is at lower energy in the compounds with Al (this band is partially filled in pure, non superconducting, $`\mathrm{AlB}_2`$). As a consequence, the extra electrons coming from Al do not fill the $`\sigma `$-bonding band completely, leaving hole pockets around $`Z`$. As we relax the structure, the $`\sigma `$ bands increase their width due to the smaller $`a`$ value, and the situation near $`E_F`$ changes only because of the occupation of the free-electron band around $`\mathrm{\Gamma }`$. These results suggest that although the rigid band scheme can be considered roughly correct, the dependence upon doping of the physical properties and of Fermi surface states relevant for superconductivity, may be crucially dependent on the actual substitutional atom.

## IV CaSi2 and related compounds

When prepared at ambient pressure $`\mathrm{CaSi}_2`$ has a rhombohedral structure (space group $`R3\overline{m}`$, $`a=3.85`$ Å, $`c=30.62`$ Å)  and is a semimetal , not superconducting down to $`30`$ mK. By annealing under pressure - typical conditions are $`8`$ GPa and $`800`$ K - a tetragonal superconducting ($`\alpha \mathrm{ThSi2}\mathrm{type}`$) phase appears with a 3-dimensional network  of Si, showing superconducivity with $`T_c=1.58`$ K. At pressures between $`7`$ GPa and $`9.5`$ GPa, rhombohedral structure samples undergo a transition to a trigonal phase ($`P3\overline{m}1`$, $`a_T3.78`$ Å, $`c_Tc_R/64.59`$ Å) with Si at $`2d`$ positions, and $`z=0.4`$, where $`z`$ is the coordinate of the B atoms in units of $`c`$ ($`z=0.5`$ corresponds to the perfect $`\mathrm{AlB}_2`$ structure). A second phase transition occurs at $`16`$ GPa. This new phase may be described as  an $`\mathrm{AlB}_2`$-like structure with $`z`$ slightly smaller than the ideal value $`1/2`$. The $`\mathrm{AlB}_2`$-type structure was previously observed in rare-earth silicides , while in CaSi2 it is apparently stabilized by doping the Si planes with Be.

In Fig. 6 we show the energy bands of CaSi2 in the $`\mathrm{AlB}_2`$ structure (dashed lines), superimposed with those of the trigonal structure (full lines), along the main symmetry lines of the Brillouin zone. Since this polymorph can be obtained only under pressure, we do not minimize the lattice parameters, but rather use the experimental values $`a=3.7077`$ Å and $`c=4.0277`$ Å, corresponding to P=19.3 GPa, and optimize the internal parameter $`z`$. The bands in the $`\mathrm{AlB}_2`$ structure show strong similarities with those of simple metal diborides, and of CaGa2. The major difference is provided by Ca $`d`$ states, which are mostly found in the conduction region, above 4 eV; however, there is a p-d mixing between Si and Ca, especially concerning Si $`\pi `$ and Ca $`d_{z^2}`$ orbitals, also mixing with the free-electron states. Unlike in MgB2, the Fermi level of CaSi2 cuts only the antibonding part of the $`\pi `$ bands, while the bonding part is fully occupied. This of course implies a reduced contribution of $`\pi `$ states to the stabilization of an $`sp^2`$ environment, consistent with the $`sp^3`$-like trigonal distortion. We will come back to this point later. The energy bands in the trigonal structure show that, a part from a small rigid shift, the $`\sigma `$ bands change very little with respect to the $`\mathrm{AlB}_2`$ structure. The $`\pi `$ bands, on the other hand, show anticrossings which in the antibonding case are located right at the Fermi level. This is most clearly indicated in the density of states of CaSi2, reported in Fig. 7. The dip at $`E_F`$ in the trigonal structure provides a textbook explanation for the lattice distortion.

In order to understand the stability of $`\mathrm{AlB}_2`$-type polymorph, we study the total energy variation as a function of $`z`$, relative to the ideal $`z=0.5`$ case, for two experimental structures, and report it in Fig. 8(a). There is a good agreement between theory and experiment (within 2.5%), and a non-negligible ($`30`$ mRy per formula unit) stabilization energy. We notice, however, that the stabilization energy decreases with pressure. As we fix $`a`$ and $`c`$ to the values corresponding to the larger pressure of 19.3 GPa, we see a larger discrepancy relative to the experimental equilibrium value of $`z`$ ($`z_{exp}0.448`$). This indicates that some different modification might be going on in the experimental samples. In order to discuss the link between the stability of CaSi2 in the $`\mathrm{AlB}_2`$ structure and the $`\pi `$band filling, without modifying the Si2 planes themselves, we vary the atomic number of the cation, which now represents a virtual atom ranging from Ca to K (corresponding to $`Z=20+x`$). Our results, reported in Fig 8(b), show that while $`x=0.3`$ leaves the total energy curve almost unchanged, $`x=0.7`$ or larger give a nearly vanishing trigonal stabilization energy. This is consistent with our previous conclusion that a filling of $`\pi `$-antibonding bands destabilizes the $`\mathrm{AlB}_2`$ structure.

This last conclusion is supported by the existence, in the hexagonal structure, of $`\mathrm{CaBe}_\mathrm{x}\mathrm{Si}_{2\mathrm{x}}`$, with $`x=0.75`$, at $`a=3.94`$ Å and $`c=4.38`$ Å. In fact, Be doping removes electrons from the Si2 planes, and brings the formal filling of the above $`\pi `$-antibonding bands to the level found in diborides. To investigate these effects, we studied the (artificially ordered, but existing) $`\mathrm{CaBe}_\mathrm{x}\mathrm{Si}_{2\mathrm{x}}`$ system for $`x=1`$, not far from the experimental stoichiometry. The corresponding bands, shown in Fig. 9, have been calculated at the experimental lattice parameters ordering Be and Si atoms as, e.g., in hexagonal BN. In this way, the Brillouin zone is the same as in the undoped compound, and the interpretation is easier. The CaBeSi band structure shown are overall similar to that of the undoped compound. The symmetry lowering induced by substitution produces the large splitting of degeneracies at the K and H points. The major difference is of course related to the position of the Fermi level, which is almost exactly coinciding with that found in the superconducting MgB2. The reduced bandwidth for the $`\sigma `$ bonding bands may be associated with the smaller size of Be orbitals. The electron depletion from the antibonding $`\pi `$-bands removes the trigonal instability, as demonstrated by total energy calculations. In fact, while the $`\mathrm{\Gamma }`$point phonons of CaSi2 in the $`\mathrm{AlB}_2`$ structure show the presence of an instable phonon, in agreement with FLAPW total energy calculations, in $`\mathrm{CaBeSi}`$ the phonons are stable, as shown in Table I. Table I shows that this system has lower phonon frequencies than the diborides, as is normal given the heavier atomic species involved. If we now immagine to vary $`\mathrm{x}`$ in $`\mathrm{CaBe}_\mathrm{x}\mathrm{Si}_{2\mathrm{x}}`$, we will therefore have that, at some value of $`x`$, the $`B_{1g}`$ mode, corresponding to the trigonal distortion, will be close to an instability. This situation may thus offer an interesting playground in the search for superconducting materials.

The Fermi level of $`\mathrm{CaBeSi}`$ cuts the $`\sigma `$ bands, leaving hole pockets quite similar to those of MgB2. While it is difficult to say, on the basis of our electronic structure calculations only, whether the similarities with MgB2 band structure can lead to superconductivity, it may be useful to provide hints on the stability of this compound, which may help experimentalists in the search for new compounds. The PDOS in Fig. 10 show a strong hybridization between Si and Be states, in the valence region. In other words, substitutions on the Si site can be possible without disrupting the $`sp^2`$ network if the one-electron energies of the impurity atoms are comparable with those of Si, and light atoms should be preferred in order to keep high values of phonon frequencies. We may suggest B, which however corresponds to a hole doping smaller than Be.

## V Conclusions

We have performed electronic structure calculations for MgB2, BeB2, $`\mathrm{Mg}_{0.5}\mathrm{Al}_{0.5}\mathrm{B}_2`$, $`\mathrm{CaSi}_2`$ in the trigonal and in the $`\mathrm{AlB}_2`$-type polymorphs at the high pressure experimental structural parameters, and for the $`\mathrm{CaBe}_\mathrm{x}\mathrm{Si}_{2\mathrm{x}}`$ system ($`x=1`$), experimentally stabilized for $`x=0.75`$. We have calculated band structures, PDOS, transport properties and phonon frequencies, using the FLAPW and the norm-conserving ab-initio pseudopotential methods, within the local density approximation. The following conclusions can be drawn on the basis of our calculations: as compared to MgB2, BeB2 has a similar filling of $`\sigma `$ bands, but differs in terms of $`\pi `$ states near $`E_F`$. The DOS at $`E_F`$ is lower, and some phonon frequencies are substantialy higher while others are comparable or lower. It would be interesting to further investigate this compound. Al substitution cannot be fully simulated by a rigid band model, as the $`\pi `$ antibonding states appear to be quite sensitive to Al substitution, and absorb much of the added electrons. Transport properties show a substantial two-dimensionality of these compounds, due to $`\pi `$ states.

We studied CaSi2 in its trigonal structure, and show that the distortion can be removed only by subtracting (in a virtual crystal-like approximation) at least $`0.7`$ electrons. In fact, $`\mathrm{CaBeSi}`$ shows no tendency towards trigonal distortion, and electron states near $`E_F`$ of similar bonding nature to those of MgB2.

## VI Acknowledgements

We thank M. Affronte and A. Gauzzi for stimulating discussions, and M.Affronte for sharing his results prior to publication. This work was partially supported by the italian Consiglio Nazionale delle Ricerche (CNR) through the “Progetto 5% Applicazioni della superconduttività ad alta Tc”.


format: text

text:
# Variable-range-hopping in two-dimensional system.

## Abstract

Computer modeling of the VRH conductivity in the two-dimensional system has been done by kinetic Monte Carlo method, which includes some new elements. Study of the temperature dependence of the conductivity, testing of the different scaling relations, and study of the size effect show that the Efros-Shklovskii mechanism of the VRH is valid in the slightest details. It has been also shown that simultaneous transitions of many electrons are not important. The reasons of disagreement with previous works are thoroughly analyzed.

The concept of the variable range hopping (VRH) of localized electrons belongs to Sir Mott. He has considered the phonon assisted tunneling of localized electrons between different sites. Mott has found that the typical hopping length for a two-dimensional (2D) case $`R_M(a/3)(T_M/T)^{1/3}`$ increases with the decreasing temperature. This is the origin of the term “VRH”. Here $`T_M=\beta _M(g_0a^2)^1`$, $`a`$ is the localization length of an electron, $`g_0`$ is the density of states at the Fermi level, $`\beta _M`$ is a numerical coefficient. The conductivity obeys the Mott’s law

$$\sigma _M(\gamma e^2/T)\mathrm{exp}\left((T_M/T)^{1/3}\right),$$
(1)
where $`\gamma `$ is a typical phonon frequency of the order of $`10^{12}\mathrm{sec}^1`$.

Electron-electron interaction has added a new chapter to the theory of VRH. The study of the interaction in localized regime has been initiated by Pollak and Srinivasan. Later on Efros and Shklovskii have shown that the single particle density of states (DS) in 2D case linearly tends to zero as energy tends to the Fermi energy. This phenomenon, which is called the Coulomb gap(CG), is due to the long range part of the Coulomb interaction which, in a sense, remains non-screened. At large disorder the DS $`G(ϵ)=2|ϵ|/\pi e^4`$ has a universal form. The DS is given by the only combination of the energy and electron charge which has a proper dimensionality. In fact, the CG results from the Coulomb law and from the discrete nature of electron charge.

Simple arguments based upon single-electron excitations show that the VRH with the Coulomb interaction (CI) obeys Efros-Shklovskii (ES) law

$$\sigma _c(\gamma e^2/T)\mathrm{exp}\left((T_0/T)^{1/2}\right),$$
(2)
where $`T_0=\beta _0e^2/a`$. The percolation approach (See references in) gives $`\beta _M=13.8,\beta _0=6.5`$. The hopping length $`R_C(a/4)(T_0/T)^{1/2}`$.

Many experimental works on the 2D structures report the ES law, but very often the data are ambiguous because one should have a really large interval of $`\sigma `$ to distinguish between the two laws. We think that the most important results have been obtained recently by groups of Jiang and Dahm and Adkins. Using an artificial screening provided by metallic electrode, parallel to a plane with the 2D electrons, they have proved that the VRH transport reflects the crucial feature of the CG, the sensitivity to the long range interaction.

The theoretical understanding of the VRH is controversial. Efros and Shklovskii argue that in 2D case only the single electron transition with the typical hopping length are important and the VRH has the same universal nature as the CG. Pollak and his followers claim that both sequential and simultaneous transitions of different groups of electrons are very important for the VRH, so that the physics of the CG which is valid for the long hops only is not relevant. These claims are based upon the theoretical idea that simultaneous transition of many electrons may provide a gain in the activation. Computer simulations of the VRH and the similar phenomena made by this group also show importance of the many-electron transport both in 3D and 2D cases but no new $`T`$-dependence of the conductivity has been proposed. The simulation of our group reveals scaling typical for the CG, but in this work a simplified model has been used which does not take into account any simultaneous transitions and underestimates sequential transitions with the short length.

We present here the results of computations which takes into account all the sequential transitions and the most important simultaneous transitions of two electrons. We show that the VRH both in the CI case and without interaction (M-case), closely obeys the laws Eqs. (2),(1) and different scaling laws which reflect the physics of the both cases. We also confirm earlier prediction that simultaneous transitions of many electrons are not important for the VRH in 2D case at low temperatures and analyze the discrepancy with the simulations by Pollak group.

We consider the standard lattice model with the Hamiltonian

$$H=\underset{i}{}\varphi _in_i+\frac{1}{2}\underset{ij}{}\frac{1}{r_{ij}}\left(n_i\nu \right)\left(n_j\nu \right),$$
(3)
where $`n_i=0,1`$ are the occupation numbers. The quenched random site energies $`\varphi _i`$ are distributed uniformly within the interval $`[1,1]`$, $`\nu `$ is the average occupation number,which is taken to be 1/2 everywhere below. The magnitude of the quenched disorder is enough to provide the universal CG at all energies important in our temperature range. Here and below the lattice constant is a length unit and the nearest Coulomb energy which is equal to the amplitude of the disorder is both energy and temperature unit. In the M-case the interaction term in Eq. (3) is neglected.

To simulate the VRH by the kinetic Monte Carlo (MC) method we add an electric field $`E`$ and consider the periodic boundary conditions (See details in Ref). To get conductivity we calculate the dipole moment due to electron transitions in the $`E`$-direction and divide it by the number of MC steps and by $`E`$. This assumes that our time unit is $`1/\gamma `$. We check that the result is $`E`$-independent. In fact we take $`E`$ from the condition $`ER_{C,M}=0.1T`$. The array size is taken $`100\times 100`$ for the temperatures above 0.05 and $`200\times 200`$ for the temperatures below 0.05. In the first approximation we neglect simultaneous transitions of many electrons. The new element of our algorithm is that we do not include tunneling probability into MC process. Instead we choose a pair of sites $`i`$ and $`j`$ with the probability $`\mathrm{exp}(2r_{ij}/a)/Z_1`$, where $`Z_1`$ is the normalization factor. Then the transition is rejected if both sites are occupied or empty, and finally it is performed with the probability $`1/(1+\mathrm{exp}(ϵ_{ij}/T))`$, where $`ϵ_{ij}`$ is the energy difference between the two configuration. This new scheme dramatically decreases the runtime of the simulation and leads to $`L^2`$, rather than $`L^4`$, algorithm in number of MC steps. Note that for the CI-case each MC step consists of $`L^2`$ operation itself due to recalculations of site energies after each transition.

To reveal and compare the physics of the VRH in both CI- and M-case we perform both simulations. Fig. 1 shows that the laws Eqs. (2),(1) are reproduced very good for both CI- and M-cases. From solid and dashed line one gets $`\beta _0=5.8`$ and $`\beta _M=13.0`$ respectively, that agrees well with the percolation values cited above. Note that the temperature where two curves become close to each other marks the upper temperature limit of the VRH. At even higher temperatures the two curves coincide and the conductivity decreases with the increasing temperature which manifests transition to a diffusion regime.

To test the origin of the VRH we check the scaling laws which reflect an intimate physics of both CI- and M-cases: $`T_0`$ and $`T_M`$ respectively are the only characteristic temperatures (energies) in the problem. These laws can be written in the universal form

$$\begin{array}{c}\sigma _M=(\gamma e^2/T_M)f_M\left(T_M/T\right),\\ \sigma _{CI}=(\gamma e^2/T_0)f_{CI}\left(T_0/T\right),\end{array}$$
(4)
where $`f_M(x)`$ and $`f_{CI}(x)`$ are some functions. It follows from the scaling laws Eq. (4) and the above definitions of $`T_M`$ and $`T_0`$ that $`\sigma _M/a^2`$ and $`\sigma _{CI}/a`$ are functions of $`Ta^2`$ and $`Ta`$, respectively. Thus, if one changes both $`T`$ and $`a`$ in the M-case in such a way that $`Ta^2`$ is constant then $`\sigma _M/a^2`$ will be also constant. For CI-case one should change $`T`$ and $`a`$ keeping $`Ta`$ constant to get $`\sigma _{CI}/a`$ constant. This test has been successfully performed and the convincing results are shown in the inset of Fig.1. Note that the small deviation of $`\sigma /a^2`$ from the constant value for the M-case is due to the lattice effect. Since $`R_M<R_C`$, the lattice effect in the M-case appears at larger $`a`$ than in the CI-case.

In Fig.2 we present the results for the size effect in the CI-case. We have calculated the average value of the conductivity av$`=\sigma _L/\sigma _{\mathrm{}}`$, where $`\sigma _L`$ is the logarithmic average of the conductivity over arrays with different disorder at given $`L`$ and $`T`$, $`\sigma _{\mathrm{}}`$ is the conductivity at the maximum $`L`$ at a given $`T`$. We have also calculated dispersion dis$`=\delta \sigma /\overline{\sigma }`$, where $`\delta \sigma `$ is the dispersion of the conductivity and $`\overline{\sigma }`$ is the arithmetic average of the conductivity at given $`L`$ and $`T`$.

We did these computations not only to get an idea how large the size of an array should be to obtain the macroscopic value of the average conductivity and what value of the dispersion one should expect at a given $`L`$ and $`T`$. Another reason is that the $`L`$\- and $`T`$-dependences of the average and the dispersion give us a new insight into the physics of the VRH. One can see from Fig. 2 that the dispersion is a function of only one argument $`TL`$, as it follows from the ES theory. Indeed, the correlation length $`L_C`$ of the percolation network is $`L_C=R_C(T_0/T)^{\nu /2}`$. Since the exponent of the correlation length $`\nu =1.33`$, one gets that $`L_C`$ is very close to $`1/T`$. Thus, $`L/L_CLT`$ which explains the $`L`$\- and $`T`$-dependence of the dispersion.

The behavior of the average is completely different. It does not change until very low values of $`TL`$. At even smaller $`TL`$ it decreases dramatically and it is not any more a function of $`TL`$. We believe that it is a function of $`L/R_{CI}L\sqrt{T}`$, and the decrease of the average conductivity appears due to the hard gap in the DS. It is known that the hard gap between occupied and empty states results from the size effect and its value is of the order of $`1/L`$. At $`L<R_C`$ the energy stripe of the VRH is completely within this gap. We suggest that this is the reason of the decrease of the average conductivity. Fig. 2 confirms qualitatively this suggestion. Unfortunately, we are unable to check it quantitatively.

Now we come to the simultaneous transitions of many electrons in the interacting system. We consider only two-electron transitions, which are the most important.

Suppose that in the initial state electrons occupy sites 1 and 3 and in the final state they are in the sites 2 and 4. The activation energy of this process $`ϵ_{1234}`$ is the energy difference between the final and the initial states. Due to Pollak it may be smaller than the sum of the energies of the single-electron transitions. Moreover, the increase in the activation exponent may beat an extra tunneling exponent, which appears because of the tunneling of two electrons. That is why Pollak and his followers insist on many-electron processes.

First we calculate the probability of the process. Let us consider the case when $`r_{12},r_{34}R`$, where $`R`$ is the distance between pairs (1,2) and (3,4). Following Ref. we concentrate on the single phonon transitions. It is obvious that these transitions only cannot provide anything but independent hops of electrons. The result is different if one considers the CI as well. We think that the amplitude of the most important process is described by the second order perturbation theory. It consists of a product of two matrix elements. The first one gives the electron transition in one of the close pairs, say (1,2), due to the phonon absorption without energy conservation. This element has a factor $`\mathrm{exp}(r_{12}/a)`$. The system comes to the intermediate state, where sites 2 and 3 are occupied. The second element describes the transition of the second electron from 3 to 4 in the external electrostatic potential created by the dipole (1,2). This element has a factor $`\mathrm{exp}(r_{34}/a)/R^3`$. The $`R`$-dependence is due to the dipole-dipole interaction. Use of the Golden rule gives the correct energy conservation law and a factor $`\mathrm{exp}(2(r_{12}+r_{34})/a)/R^6`$. The result is the same as for the Förster diffusion of molecular excitons , very similar but not identical problem. The transition rate used in Ref. and many others is taken to be $`R`$-independent at large $`R`$. We cannot understand this result because to get the common energy conservation law one should assume that the remote pairs can “talk” to each other, which means there should be an interaction between them with some $`R`$-dependence. We argue below that the contradiction between our results originates from the $`R`$-dependence.

To accommodate this fast decay with $`R`$ we assume that the simultaneous transitions occurs only if the initial or final position of the second electron is within the spheres with the radii $`R_C`$ from the initial or final position of the first electron. Note that this assumption overestimates the role of the double electron transitions in the most important case when one electron makes a long hop $`R_C`$ and the other one makes a short hop. The interaction between such pairs decreases (as $`1/R^4`$) at distances smaller than $`R_C`$.

To simulate the double hop transport we use the same idea as for single hop case. This time, however, we choose a pair of sites $`i`$ and $`j`$ for single hop transition with the probability $`p_1=\mathrm{exp}(2r_{ij}/a)/Z_2`$, and sites $`i`$, $`j`$, $`k`$ and $`l`$ for double hop with the probability $`p_2=\mathrm{exp}(2(r_{ij}+r_{kl})/a)/Z_2`$, where $`Z_2`$ is found from the normalization condition for the sum of two probabilities.

If the double hop is chosen, we check the occupation for both pairs of sites $`(i,j)`$ and $`(k,l)`$ and then use the same Boltzman factor as before with $`ϵ_{ijkl}`$ instead of $`ϵ_{ij}`$. The algorithm is also $`L^2`$ in the number of MC steps. However, since the number of possible transitions of two electrons is much larger than the number of single transitions, one can show that there are $`4R_C^2\mathrm{exp}(2/a)`$ attempts of double transition per one attempt of a single one, which makes simulation much more time consuming.

The contribution of single and double hops to the total conductivity can be separated by collecting their dipole moments in the different cells. We have checked that the single hop contribution does not change in the presence of double hops. The results of the simulation are shown on the inset of Fig.2. We are able to calculate double hop contribution for the case $`a=0.5`$, because at larger $`a`$ the value of $`R_C`$ is too large. The single hop contributions are shown for both $`a=1,0.5`$. Note that a small deviation from the universal behavior at high temperatures is a result of the lattice effect for $`a=0.5`$. One can see that the double hop contribution is from one to two orders of magnitude smaller than the single hop contribution at VRH region of temperatures. Moreover, the relative fraction of double hops in the total conductivity decreases with decreasing temperature. This result does not support Pollak’s prediction of the gain in the activation energy for simultaneous transition of two electrons as compared with the single electron transition.

The computational data, obtained in the papers  and many other works of this group show the importance of the many-electron hops which contradicts to our result. In our opinion the $`R`$-independent probability of the simultaneous hops assumed by this group is responsible for the contradiction.

Indeed, if in a large array each pair has an equal probability to make a transition with any other $`L^2`$ soft pairs in the system the number of transitions for the pair per time is proportional to $`L^2`$. This is a giant overestimate of the role of many-electron transitions. The abnormal size effect should be a manifestation of this mistake. Unfortunately, we have found only one work, where authors take care about the size effect and this work deals with the 3D case. However, we think that the infinite range interaction between pairs should provide the same effect in the 3D-case. Looking at Fig. 3 of Ref. one can see an exponential size effect without any signs of a saturation. The resistivity increases with the size because the algorithm include the selection of states with the highest transition rate among them. The selected states must consist of soft pairs which are far from each other and do not contribute to the dc conductivity.

Another reason for contradictions is that the method used in all these papers is applicable to the mesoscopic systems only. It consists of selecting a number of low-energy states of the total system to study transitions between them. It looks nice because it takes into account all many electron transitions, while we could afford the simultaneous transitions of two electrons only. The problem is in the necessary number of states, which provides reasonable thermodynamic and kinetic description. This number increases exponentially with the number of sites in the array. Due to the sensitivity of the CG to the long range interaction the size of the array $`L`$, which provides macroscopic regime, increases with decreasing $`T`$. One can estimate from Fig. 2 that $`L3/T`$ gives a reasonable dispersion. For example, T=1/400 is the middle of the temperature interval in Ref.. It follows that $`L`$ should be at least 1000. Instead, authors compute few samples with different disorder at $`L22`$. These samples have very different values of conductivities and different patterns of the $`T`$-dependences (See Fig.1 of Ref.). This is a typical mesoscopic result.

Finally, our computer simulation confirms the ES theory of the VRH transport in 2D systems and shows that the attempts of its revision have no ground.

We are grateful to O. Biham and B. I. Shklovskii for the interest to this work and to E. I. Rashba for thorough discussion of the theory of simultaneous transitions. The work have been funded by the US-Israel Binational Science Foundation Grant 9800097, the computations have been made in CHPC of the University of Utah.


text:
# Dynamic Response of Adhesion Complexes: Beyond the Single-Path Picture

## Abstract

We analyze the response of molecular adhesion complexes to increasing pulling forces (dynamic force spectroscopy) when dissociation can occur along either one of two alternative trajectories in the underlying multidimensional energy landscape. A great diversity of behaviors (e.g. non-monotonicity) is found for the unbinding force and time as a function of the rate at which the pulling force is increased. We highlight an intrinsic difficulty in unambiguously determining the features of the energy landscape from single-molecule pulling experiments. We also suggest a class of “harpoon” stickers that bind easily but resist strong pulling efficiently.

The last decades have witnessed a remarkable development of physical investigation methods to probe single molecules or complexes by various micromanipulation means. New techniques have been put forward to probe the unfolding of proteins and to quantify the strength of adhesion structures . An important step in this direction is the proposal of the group of Evans to use soft structures to pull on adhesion complexes or molecules at various loading rates (dynamic force spectroscopy) . Moving the other end of the soft structure at constant velocity induces on the complex a pulling force that increases linearly in time $`f=rt`$. Measuring the typical rupture time $`t_{\mathrm{typ}}`$ yields a typical rupture force $`f_{\mathrm{typ}}=rt_{\mathrm{typ}}`$ that depends on the pulling rate $`r`$. This provides information as to the energy landscape of the bound complex. Indeed, in many situations one observes a linear increase of $`f_{\mathrm{typ}}`$ with $`\mathrm{log}(r)`$, which can be understood within a simple adiabatic Kramers picture for the escape from a well (bound/attached state) over a barrier of height $`E`$ located at a projected distance $`x`$ from the well along the pulling direction. The progressive increase of the force results in a corresponding increase of the escape rate, so that, in agreement with some experiments , the typical rupture force increases logarithmically with $`r`$: $`f_{\mathrm{typ}}k_\mathrm{B}T/x\mathrm{ln}[rx/(k_\mathrm{B}T\omega )]`$, where $`\omega `$ is the escape rate in the absence of force. The rupture time on the other hand decreases with $`r`$. The occurrence in some cases of two successive straight lines in a $`[f_{\mathrm{typ}},\mathrm{log}(r)]`$ plot has been argued to be the consequence of having two successive barriers along the 1D escape path, the intermediate one showing up in the response at fast pulling rates  (Figs. 1a and 2). Other theories have tried to back up more complete information as to the overall effective 1D potential landscape by an analysis of the probability distribution for rupture time and of the statistics of trajectories before rupture . Assemblies in series and in parallel of such 1D bonds have also been considered .

In this Letter we point out limitations arising from the a priori assumption of a single-path topology of the energy landscape for the interpretation of such experiments. From the analysis of simple examples with a two-path topology, we draw three conclusions: (i) first, the dependence of the rupture force and rupture time on the pulling rate can take various forms, including non-monotonic behavior (see e.g. Figs. 3 to 5). (ii) Second, the main features of the energy landscape can not be unambiguously deduced from a \[$`f_{\mathrm{typ}},\mathrm{log}(r)`$\] plot, as very different landscapes can yield similar curves (Fig. 6). (iii) Third, we propose simple “harpoon” designs (Fig. 1 c and d) for functionally efficient stickers that can bind easily but resist strongly in a range of pulling forces (Fig. 4).

Obviously for real binding/adhesion complexes, there are numerous (conformational) degrees of freedom, and the configurational space is clearly multidimensional. This allows for complex energy landscapes and various topologies for the structure of their valleys and passes . Only the probing (pulling) is unidirectional. We note in passing that even for more macroscopic sticky systems, usual adhesion tests for soft adhesives often show up hysteresis loops associated with the existence of more than one degree of freedom . We do not attempt here an exhaustive exploration of effects allowed by the multidimensionality of the phase space, but rather focus on a few simple two-path topologies (Fig. 1), to argue for the three points mentioned above.

The three examples we consider, sketched in Figure 1 b, c, and d, correspond to simple hairpin schemes whereby detachment can proceed through two alternative routes $`\alpha `$ and $`\beta `$. These simple quasi 1D schematic situations can be conveniently dealt with using an adiabatic Kramers theory, which has been shown to be an efficient way of obtaining semi-quantitatively correct answers .

A common set of notations can be ascribed for all cases (Fig. 1). From the fundamental bound state “0”, the route $`\alpha `$ for escape (detachment) is over barriers $`a`$, of height $`E_a`$ located at a projected distance $`x_a`$ from “0”. Alternatively, escape can occur through branch $`\beta `$, over barrier $`b`$, of height $`E_b`$ and projected distance $`x_b`$. All energies and projected distances are measured relative to the state “0” (i.e. $`E_0=0`$ and $`x_0=0`$). Intermediate barriers $`a^{}`$, $`b^{}`$ and local minima A and B may exist, with energies $`E_a^{},E_b^{},E_A,E_B`$ (all positive), and projected distances $`x_a^{},x_b^{},x_A,x_B`$. In line with typical values from experiments, we choose to write energies in units of $`k_\mathrm{B}T4`$ pN nm and distances in nm.

Practically, we describe the time evolution of the probabilities of being in the potential minima (bound states) using “chemical” transition rates over the barriers as given by the Kramers formula. We furthermore assume the attempt frequencies to be constant and all equal to $`\omega _0`$ which provides the only intrinsic time-scale in the problem, so that the transition rate from minimum $`I`$ over the neighboring barrier $`i`$ is $`\omega _0\mathrm{exp}[(E_iE_I)+f(t)(x_ix_I)]`$. For the plots of Figures 2 to 6 we take arbitrarily $`\omega _0=10^8`$ s<sup>-1</sup>. Jump over the rightmost barrier ($`a`$ or $`b`$) of either path corresponds to rupture leading to escape to $`x\mathrm{}`$.

We focus on the case where either $`E_b^{}`$ or $`E_b`$ is larger than $`E_a`$ so that $`\alpha `$ is the “natural” route by which attachment and detachment proceeds in the absence of pulling. We also limit ourselves to simple scenarios in which the force is linearly increased in time $`f=rt`$.

For further reference we recall the classical single-path scenario (Fig. 1a) in Figure 2 for a typical set of parameters, and then we turn to a brief analysis of the three geometries we have introduced (Fig. 1 b, c, and d).

First case: switch – Topology as in Figure 1b. Escape occurs through either barrier $`a`$ or barrier $`b`$ both located downwards in the pulling direction ($`x_a,x_b>0`$). The escape proceeds through path $`\alpha `$ at weak pulling rates as $`E_a<E_b`$, but if $`x_b>x_a`$ it can switch to path $`\beta `$ for pulling forces $`f`$ large enough such that $`E_afx_a>E_bfx_b`$. The result (see Fig. 3) is then a succession of two straight lines of decreasing slopes in the \[$`f_{\mathrm{typ}},\mathrm{log}(r)`$\] plot, the first one (slope $`1/x_a`$) characteristic of the spontaneous route $`\alpha `$ while the second (slope $`1/x_b`$) provides information on the alternative route $`\beta `$. In the trivial case $`x_a>x_b`$ route $`\beta `$ is never explored so that the classical single-path picture applies.

To clarify the calculation leading to the plot in Figure 3, we describe the evolution of the probability of attachment $`p(t)`$ at time $`t`$ of the system initially attached at time $`t=0`$ \[$`p(0)=1`$\] by

$$_tp(t)=\omega _0(e^{E_a+f(t)x_a}+e^{E_b+f(t)x_b})p(t)$$
(1)
Solving (1) numerically with $`f=rt`$ yields $`p(t)`$ and therefore the probability density $`P(f)=\frac{1}{r}_tp(\frac{f}{r})`$ for the unbinding force. The typical values of $`f`$ are highlighted in the plots, with the whole distribution $`P(f)`$ suggested through a grey-scale. Similar procedures will be used in the following examples, with thermal equilibrium between the bound states assumed as initial conditions.

Second case: harpoon – Topology similar to the previous one but with $`x_a<0`$ (Fig. 1c). The main feature here is that as the pulling force increases, the probability to escape over $`a`$ decreases. Therefore the system gets “stuck” in route $`\beta `$. If the barrier $`E_b`$ is infinite (left side of Fig. 1c), there is a finite probability $`p_{\mathrm{}}=\mathrm{exp}(\frac{\omega _0e^{E_a}}{r|x_a|})`$ that unbinding never occurs. For a finite but high barrier $`E_b`$, pulling eventually results in unbinding but at high rupture forces (see Fig. 4). The topology thus allows here to form “easily” (i.e. over barrier $`a`$) a “harpoon” sticker that can resist strong pulling. Correspondingly the mean unbinding time increases first with pulling rate (a phenomenology connected to the negative resistance analyzed in Ref. ), before decreasing for larger values when activated escape over $`b`$ dominates. Note that the probability distribution $`P(f)`$, now consists of two separate ensembles, which coexist over a narrow region of pulling rates. This is in contrast with Figure 3 where there is a continuous evolution of a single cloud.

Third case: combo – The alternative route consists of two barriers and a local minimum B (Fig. 1d), and we focus on the case where $`E_b^{}`$ is smaller than the two others.

Thanks to the increased complexity and number of parameters in this case many scenarios can occur, covering features already unveiled in Figures 2 to 4 (e.g. switch and harpoon). More intricate pictures can also show up, as depicted in Figure 5. An explanation of this example is given in the caption, illuminating how for low or high pulling rates barrier $`a`$ controls the behavior, whereas for intermediate values, the secondary and stronger barrier $`b`$ limits unbinding. Two features are striking. First, the unbinding force (typical or average) is no more monotonic. Second, branch $`\beta `$ results in a strengthening of the binding complex for a given window of pulling rates $`r`$ (selective harpoon).

Discussion – With the three simple examples above, we have clearly enlarged the numbers of behaviors that one may obtain from a classical dynamic force spectroscopy method (see Figs. 3 to 5). Conversely, we also want to stress the second point (ii) mentioned in the introduction: simple patterns (e.g. the succession of two lines of increasing slopes) can be the outcome of many diverse landscapes. For example, Figure 6 displays force-rate curves similar to that of Figure 2, but that correspond to sensibly different landscapes. Not only are the typical and average unbinding forces very similar in the three cases, but so are the probability distributions for most values of $`r`$. Only close to the cross-over between the two straight lines can slight differences be detected. To distinguish more selectively possible landscapes, it may be necessary to use other temporal sequences than the simple $`f=rt`$, e.g. to reveal equilibration processes between local minima.

Eventually we would like to emphasize that the harpoon geometries proposed here constitute a very obvious paradigm for efficient stickers. Attachment of the sticker can proceed through route $`\alpha `$ with a possibly not too high barrier $`E_a`$. The harpoon configuration then allows to benefit from the much stronger $`b`$ barrier for a given window of pulling forces, making the sticker more efficient in these conditions. This “hook” design is obviously a favorable strategy for adhesion complexes, the function of which is to maintain adhesion under the action of well-defined tearing stresses. It would be surprising if advantage was not taken of this by some biological systems.


text:
# A microscopic model for exchange bias

## Abstract

The domain state model for exchange bias has been further investigated by considering vector spins for the antiferromagnet instead of Ising spins used in the earlier studies. The qualitative results are similar to those with infinite anisotropy for the antiferromagnet. However, under certain conditions softer spins can lead to an even stronger bias field. The study shows a nontrivial dependence of the exchange bias on the antiferromagnetic anisotropy.

Theoretische Tieftemperaturphysik, Gerhard-Mercator-Universität Duisburg, 47048 Duisburg, Germany.

Keywords: Exchange biasing, Domain structure, Dilute antiferromagnet, Heisenberg antiferromagnet

The shift of the hysteresis loop along the field axis in a system containing a ferromagnet(FM)-antiferromagnet(AFM) interface is known as exchange bias (EB). The discovery of this effect by Meiklejohn and Bean in Co-CoO system  dates back to 1956. Subsequently the effect was observed in a variety of different materials (see  for a review). There have been several approaches towards a theoretical understanding of the phenomenon and numerous different mechanisms are believed to be responsible for the shift. Malozemoff  explained the shift as a result of interface roughness. Koon considered a spin-flop coupling between the FM and the compensated AFM interface as responsible for the effect . However, Schulthess and Butler argued that the spin-flop coupling alone cannot account for the effect  but an interface coupling through uncompensated defects together with the perpendicular coupling could explain the mechanism .

Recently a model was proposed  where the cause of the shift was attributed to the extra exchange field provided by the remanent magnetization of the domains in a diluted AFM. The model explained successfully several aspects related to the effect. The essential idea behind the model is that when the diluted AFM is cooled below the N$`\stackrel{´}{\text{e}}`$el temperature in presence of a field, it ends up in a domain state as opposed to long range order. These domains are metastable since the domain walls require large thermal fluctuations for any subsequent movement. The AFM is diluted by replacing its atoms by nonmagnetic impurities (e.g. $`\text{Co}_{1x}\text{Mg}_x\text{O}`$). The domains are usually pinned at the impurity sites because at those sites it does not cost any energy for the creation of a domain wall. These domains carry a remanent magnetization  which provides the biasing field to the FM causing the shift of the hysteresis loop. Several issues related to EB, such as temperature dependence, training effect , role of AFM thickness  and, of course, the dependence of dilution on EB have been successfully explained within the framework of this model. In the earlier work , owing to the strong uniaxial anisotropy of CoO, the AFM was described by an Ising model. In the present work, we relax this restriction on the AFM. In order to investigate a broad class of systems for the AFM a Heisenberg model with variable uniaxial anisotropy is considered.

The system we will study consists of one monolayer of FM and $`t`$ monolayers of AFM. A certain fraction $`p`$ of the AFM sites are left without a spin. The FM is exchange coupled to the topmost layer of the AFM. The Hamiltonian is thus given by

$``$ $`=`$ $`J_{\mathrm{FM}}{\displaystyle \underset{i,j}{}}\stackrel{}{S_i}\stackrel{}{S_j}{\displaystyle \underset{i}{}}\left(d_zS_{iz}^2+d_xS_{ix}^2+\stackrel{}{S_i}\stackrel{}{B}\right)`$
$`J_{\mathrm{AFM}}{\displaystyle \underset{i,j}{}}ϵ_iϵ_j\stackrel{}{\sigma _i}\stackrel{}{\sigma _j}{\displaystyle \underset{i}{}}ϵ_i\left(k_z\sigma _{iz}^2+\stackrel{}{\sigma _i}\stackrel{}{B}\right)`$
$`J_{\mathrm{INT}}{\displaystyle \underset{i,j}{}}ϵ_j\stackrel{}{S_i}\stackrel{}{\sigma _j}`$
where $`\stackrel{}{S_i}`$ and $`\stackrel{}{\sigma _i}`$ denote spins at the $`i`$th site corresponding to FM and AFM respectively. The first line of the Hamiltonian describes the contribution to the energy coming from the FM with the $`z`$-axis as its easy axis $`\left(d_z>0\right)`$ and the $`x`$-axis as its hard axis $`\left(d_x<0\right)`$. The in-plane anisotropy keeps the FM preferentially in the $`yz`$ plane. The second line is the contribution from the AFM with quenched disorder $`\left(ϵ_i=0,\mathrm{\hspace{0.25em}1}\right)`$ also having its easy axis along $`z`$ $`\left(k_z>0\right)`$. The last term describes the interaction of the FM with the interface AFM monolayer. The whole system is placed in an external magnetic field $`\stackrel{}{B}=B\widehat{z}`$. In our simulations we set $`J_{\mathrm{FM}}=2J_{\mathrm{AFM}}=2J_{\mathrm{INT}}`$.

We use Monte Carlo methods with a heat bath algorithm. Each spin is at every Monte Carlo step subject to a trial step consisting of a small deviation from the original direction followed by a total flip. This two-fold trial step can take care of a broad range of anisotropies starting from very soft spins up to the Ising limit . The lateral dimension of the system was chosen in such a way that multiple AFM domains can be observed which is the case for a lattice size $`128\times 128\times (t+1)`$. Starting from a temperature above the N$`\stackrel{´}{\text{e}}`$el temperature $`\left(T_N\right)`$ of the diluted AFM but below the Curie temperature $`\left(T_c\right)`$ of the FM, the system is cooled below $`T_N`$ in presence of an external magnetic field $`\stackrel{}{B}=0.25J_{\mathrm{FM}}\widehat{z}`$. Then the hysteresis curve of the system is calculated along the $`\widehat{z}`$ direction. The EB field is determined as $`B_{\mathrm{EB}}=\left(B^++B^{}\right)/2,`$ where $`B^+`$ and $`B^{}`$ are the values of the field where the $`z`$-component of the magnetization of the FM vanishes when the field is decreasing and increasing, respectively.

From the simulations the following qualitative results emerge. In absence of any anisotropy in the FM we observe a perpendicular coupling between FM and AFM, only at very low dilution of the latter. The magnetization reversal in the FM is by coherent rotation. The picture changes with increasing uniaxial anisotropy in the FM and upon further dilution of the AFM. The magnetization reversal in the FM is now by domain wall motion and the perpendicular coupling becomes less significant. This is because uniaxial anisotropies in both the FM and the AFM having the same axis no longer lead to an energy minimum at perpendicular coupling across the interface. Moreover AFM spins with missing AFM neighbours can lower their energy by rotating parallel to the its FM neighbour. Therefore, in the framework of our calculations a spin-flop coupling is not an important mechanism for the shift of the hysteresis loop.

Coloured pictures of typical staggered domain configurations of the bulk AFM with different anisotropies can be found in the following link: http://www.thp.uni-duisburg.de/~arko/picture.html. The snapshots are taken when the net magnetization of the FM is nearly zero during the decreasing field part of the hysteresis loop. As expected the domain walls are wider for smaller values of $`k_z`$, whereas Ising like domain structures are obtained for stronger anisotropies. Typical hysteresis loops for the FM and the interface monolayer of the AFM are shown in Fig. 1. The upward shift of the hysteresis loop for the interface layer of the AFM proves the existence of remanent magnetization in the AFM domains. The magnitude of the EB field strongly depends on the amount of this upward shift.

We have calculated the EB field for a wide range of values of $`k_z`$, starting from very soft spins to rigid, Ising-like spins. Fig. 2 shows result for different thicknesses of the AFM for two values of the dilution, $`p=0.4`$ and $`0.6`$. A logarithmic scale is used along the $`k_z`$ axis only to increase the clarity of the figure. Fig. 2(a) shows that for $`p=0.4`$ a thick AFM produces a peak in the EB field at an intermediate value of $`k_z`$ while at lower thicknesses of the AFM the EB field increases with the anisotropy and saturates in the Ising limit. Qualitatively this can be understood as follows. The AFM domains are required to carry surplus magnetization at the interface along the $`\widehat{z}`$ direction which must be stable during the hysteresis loop in order to produce any EB. At low values of $`k_z`$ the minimum energy configuration of the AFM spins is predominantly dictated by the FM and the external field. In this limit the AFM will follow the FM, thereby producing no bias field. By increasing the anisotropy the AFM spins are forced to lie more along the $`z`$-axis and hence the bias increases. However, there exists a counter effect. Imagine a domain is formed upon field cooling which preferentially goes through defects thereby minimizing the energy. If now the anisotropy $`k_z`$ is increased the domain wall width will decrease and the energy to create such a wall will increase. Thus the system will respond by flattening the domain boundaries. This reduces the probability for the domains to carry any surplus magnetization and hence the bias decreases for larger values of $`k_z`$. The compromise between these two opposite effects is achieved at some intermediate value of $`k_z`$ where the bias shows a peak. However, the peak disappears at lower values of $`t`$ (fig. 2(a)) or at a higher dilution (fig. 2(b)). This happens when we are close to the percolation threshold where the domain walls pass nearly exclusively through the defects costing very little or no energy. Since no smoothening of the domain walls is required now to minimize the energy of domain formation, the second mechanism discussed above does not come into play. Hence exchange bias increases with $`k_z`$ till it saturates in the Ising limit. Since the percolation threshold in three dimensions is higher than that in two dimensions, increasing the number of AFM monolayers at $`p=0.4`$ we go below the threshold and the peak in the bias field reappears.

In conclusion, we find that the domain state model for EB proposed originally for the Ising AFM is not restricted to this limit. Rather under certain combination of thickness and dilution of the AFM, the softness of the AFM spins can lead to an even stronger bias field. There are several mechanisms responsible for the bias field, such as dilution, thickness and anisotropy of the AFM. Although a qualitative understanding regarding the dependence of EB on these parameters has been achieved, a quantitative study of the domain structure both at the interface and in the bulk of the AFM would provide a deeper understanding to the problem.

This work has been supported by the Deutsche Forschungsgemeinschaft through SFB 491 and Graduiertenkolleg 277.

Figure Captions

Fig. 1. Typical hysteresis loop along $`\widehat{z}`$ of (a) the FM and (b) the interface monolayer of AFM, for $`p=0.6`$ and $`t=2`$. The net magnetization is shown in units of the saturation magnetization and the field in units of $`J_{\mathrm{FM}}`$.

Fig. 2. Dependence of the exchange bias field, shown in units of $`J_{\mathrm{FM}}`$, on the AFM anisotropy for (a) $`p=0.4`$ and (b) $`p=0.6`$ at different thicknesses.

Figure 1: Misra et al.

Figure 2: Misra et al.


text:
# Instability of scale-free networks under node-breaking avalanches

## Abstract

The instability introduced in a large scale-free network by the triggering of node-breaking avalanches is analyzed using the fiber-bundle model as conceptual framework. We found, by measuring the size of the giant component, the avalanche size distribution and other quantities, the existence of an abrupt transition. This test of strength for complex networks like Internet is more stringent than others recently considered like the random removal of nodes, analyzed within the framework of percolation theory. Finally, we discuss the possible implications of our results and their relevance in forecasting cascading failures in scale-free networks.

The recent months have witnessed a great effort devoted to the unraveling of the properties of complex networks . These properties include the rules followed in the process of formation of the net, the resulting connectivity distribution of the networks and their robustness under unfavorable circumstances like the presence of acting damaging agents, etc . This general interest is mainly due to the fact that the subject of complex networks has a considerable impact on many branches of science and technology and also in sociology . An important observation has been to recognize that some significant networks have a scale-free connectivity distribution which means that the number of links emerging from one node statistically follows a power-law distribution $`P_kk^\gamma `$. In particular, it has also been noted that the Internet belongs to this class of networks with several studies reporting an exponent $`\gamma =2.2\pm 0.1`$ .

In this Letter, we explore the robustness of large scale-free networks in a scenario in which the failure of a node may trigger the subsequent failure of its neighbors. Two nodes are considered as nearest neighbors if they are connected by a direct link. It is thus clear that the possibility of having $`2^{nd}`$-step failures can in turn induce $`3^{rd}`$-step failures, etc and thus a breaking avalanche can be generated. This idea is not just an entertainment of theoreticians; as is now well-known, on 10 August 1996, a fault in two power lines in Oregon led, through an avalanche of cascades, to a large-scale blackout in the U.S. and Canada.

To implement our model, we will use the framework of the so-called fiber-bundle models (FBM) . In FBM a set of $`N1`$ fibers (elements) is located on the sites of a -usually regular- lattice, and one assigns to each element a random strength threshold sampled from a given probability distribution (the Weibull distribution is frequently used for this purpose). Then, the set is loaded by uniformly elevating the weight acting on each element up to a certain value, $`\sigma `$. All the elements whose thresholds are lower than $`\sigma `$ fail in the first instance. The individual load carried by each of the broken elements is equally transferred to their surviving nearest neighbors and therefore the rupture of an element may induce secondary failures which in turn may trigger tertiary failures and so on. These systems are usually conservative which means that when equilibrium is finally attained, i.e., when there occur no more casualties, the total weight borne by the $`N_s`$ surviving elements is equal to $`N_s\sigma `$. The translation of the FBM terminology to that of complex networks studied here is simple. A fiber may be viewed as a node of the net, the directions of the load transfers are now the links of the net connecting the nodes to each other and the load can represent, for instance, the intensity of electric current flowing into the nodes of the network, or the viral pressure in Internet, etc, i.e. any magnitude or agent able to surpass a security threshold of an individual node, break it, and then provoke an increase of potential damage in the neighborhood of that node.

Recently, several scenarios of instability have been considered for complex networks. In particular for scale-free networks the random removal of nodes in a fixed proportion and its impact on the global connectivity and functionality of the network has been explored  with tranquilizing conclusions: networks with $`\gamma <3`$ are completely stable. In this scenario, akin to percolation, the removal of a node never leads to cascades as it does in our model. This implies that the test of strength analyzed here is more stringent and thus offers more guarantees for the security of networks with power law connectivity distributions.

To study the instability of scale-free networks under node-breaking avalanches, we first construct a network using the Barabási-Albert (BA) algorithm . This is a stochastic growth model in which, at each time step, a new node is added and attached preferentially to the already existing ones. At the initial state, we start from a small number $`m_0`$ of disconnected nodes and the network grows by adding one new node at each time step. This new node is connected preferentially to $`m`$ old nodes with a probability that depends on the node connectivities $`k_i`$ through the relation $`\mathrm{\Pi }(k_i)=k_i/_jk_j`$. By iterating this scheme a sufficient number of times, a network consisting of $`N`$ nodes with connectivity distribution $`P_k=2m^2k^3`$ and average connectivity $`k=2m`$ develops . It is worth noting that the preferential attachment rule introduced in the BA model accounts for the rich-gets-richer feature observed in many real complex networks. Besides, the BA model has been recently improved by adding several new ingredients in order to account for connectivity distributions with exponents $`2<\gamma <3`$. The main feature of the scale-free networks is that each node has a statistically significant probability of having a very large number of connections compared to the average connectivity of the network $`k`$. This is not the case for other complex networks  where the connectivity distribution is peaked at $`k`$ and decays exponentially fast for $`kk`$ and $`kk`$. Thus we expect that working with $`\gamma =3`$ does not alter the results one would obtain for the general case $`2<\gamma 3`$, a fact confirmed by preliminary numerical simulations in more general SF networks .

Let us now assume that the scale-free network previously generated is exposed to an external pressure or force $`F`$ and that each node of the network represents an individual element able to support a finite amount of “load” $`\sigma `$. As noted before this could be seen as a system where the individual elements are continuously subjected to external agents able to affect their functionality if they overcome a given security threshold. We will also assume that in the initial state this external force is equally distributed among all the nodes in the network so that each element bears a load $`\sigma =F/N`$. Furthermore, we assign to each node a statistically distributed security threshold $`\sigma _{i_{th}}`$ ($`1iN`$) taken from a probability distribution. If the load acting on a node surpasses its threshold, the node fails and its load is equally transferred to the non-failed nodes directly linked to it. This may provoke other nodes to collapse and the cascades of failure events last until all the sound elements in the network bear a load lower than their threshold values. In order to assign the threshold values we will use the Weibull distribution $`P(\sigma _{i_{th}})=1e^{(\sigma _{i_{th}})^\rho }`$, where $`\rho `$ is the so-called Weibull index which controls the degree of threshold disorder in the system (the bigger the value of $`\rho `$, the narrower the range of threshold values). This allows us to compare the stability of systems having different levels of heterogeneity in their security threshold distribution.

The threshold rule introduced above has been used for many years to study a wide class of nonequilibrium phenomena . However, they have been extensively studied mainly for regular lattices. In the present model, the cascading process not only depends on the thresholds of the elements but also on the distribution of the neighbors, i.e., a casualty at a node is determined by both the threshold of that node and the number (and state) of nodes directly linked to it.

We have performed large scale numerical simulations of the cascading process produced by repeatedly applying the rules stated above. In the initial state all the elements that form the network are subjected to a small individual load $`\sigma `$. If that load is bigger than one (or several) threshold(s), a cascading process could start that lasts until the system arrives at a new equilibrium state where all the nodes support a load lower than their security thresholds. Then, the complete process is repeated again by imposing in the initial state (with all the elements in the non-failed state) a bigger load and applying the same failure rules. After each cascading event, the damage to the system increases, which affects both the properties of the network and its functionality. Each simulation is performed many times to average over the security threshold distribution and in the end a kind of phase diagram for each magnitude characterizing the final damaged state of the system is obtained.

The results obtained indicate that the system has an abrupt transition in its connectivity. Because of the similitude between this behavior and the one found in fracture systems on regular networks, we would call it a critical point. Figure 1 shows the size of the giant component in a BA network composed by $`N=10^5`$ nodes as a function of the control parameter $`\sigma `$ for a Weibull distribution with two different values of $`\rho `$ and for a uniform distribution of thresholds. The size of the giant component, which plays the role of an order parameter, is defined as the total number of intact nodes remaining in the largest component of the network after the cascading events divided by the system size $`N`$. We performed a depth-first search for the largest component  and averaged the results after many realizations to get a stable mean value of this component’s size. As can be seen from the figure, for small values of the weight imposed over the system, the network remains almost intact and the giant component size is still large enough to ensure the network’s functionality. However, as the load is increased, the cascading failure begins to reach more and more nodes up to the critical point where the size of the giant component suddenly falls close to zero provoking the rupture of the system in many small parts losing its properties as a functional network.

Figure 1 also shows that the network with the highest degree of homogeneity in the threshold distribution is more resilient to breakdown. This behavior is the opposite to what is seen in regular networks. Although the critical load at which the network loses its functionality shifts to the right as the level of homogeneity in the thresholds is increased, the precursory activity is less intensive and so the final breakdown of the network arises more abruptly and catastrophically, without previous warning. This result agrees with what is seen in fracture processes. Obviously, from a practical point of view, this is as unwanted as having a low critical value that makes the network very unstable. An intermediate value could satisfy both criteria; one is to have a robust network and the other to guarantee that the failure of the network is preceded by an important precursory activity which helps to foresee the cascades and the imminent collapse.

Another way to shed light on the cascading process is to inspect the change in the topology of the network as the control parameter varies. The transition from a functional network to the fragmented one is illustrated by the simulation results shown in Fig. 2, where we plot the probability that a node has connectivity $`k`$ when the system has reached the final equilibrium state. For very low damage pressures (for instance, $`\sigma =0.05`$ in the figure), the topology of the network remains unchanged. Right after the critical point ($`\sigma =0.52`$ in Fig. 2), the system loses its properties as an scale-free network and the whole system becomes disconnected, the largest connectivity of the nodes being of about 4 links.

We have also monitored the nodes that collapse more frequently. This can be observed in Fig. 3 where we have represented the fraction of broken nodes $`n_k`$ with respect to their connectivities for the same parameters of Fig. 2. The existence of a critical point is again clearly demonstrated. As the load supported by the system approaches the breakdown point from below, the nodes with higher connectivities are the most affected by the cascading process although the system still conserves a few hubs that make possible the existence of a large giant component that keeps a significant fraction of nodes interconnected. Right after the critical point only nodes forming isolated clusters remain, none of them having a large connectivity. While the fragility of scale-free networks with respect to the removal of highly connected nodes has been recently documented by several authors , it is worth noting that the above result, although pointing in the same direction, is obtained as a consequence of the model’s rules instead of being imposed from the outside, i.e., the system is subjected to an external force and it evolves according to a simple threshold rule in contrast to other models where the removal of hubs is directed.

Finally, we have characterized the cascading process itself by measuring the avalanche distribution. The size $`s`$ of an avalanche is defined as the total number of nodes that break simultaneously. The cumulative distributions $`P(s)`$ of avalanche sizes for a network composed of $`10^5`$ nodes have been represented in Fig. 4 for two different values of the threshold disorder parameter. In both cases, the avalanche size distribution was measured at the critical points which are in these cases $`\sigma _c=0.52`$ and $`\sigma _c=0.37`$ for $`\rho =5`$ and $`\rho =2`$ respectively. The distributions can be fitted, for low values of $`\rho `$, by a power law of the form $`P(s)s^\tau `$, giving an exponent for the probability distributions of $`\tau +1=1.12\pm 0.03`$. As the threshold distribution gets more homogeneous (i.e., bigger $`\rho `$), the avalanche size distribution becomes a strongly peaked function around a large mean value. This occurs because, as we move to the region of large critical points values (increasing $`\rho `$), there is a very poor precursory activity and almost all the nodes forming the network breakdown in only one time step. Above the critical point, the probability distribution of avalanches splits into two parts and is characterized by the excess of large avalanches signaling that we are in the supercritical region. Without providing a figure, we also want to report that in scale-free networks the value of the critical stress, for the same $`\rho `$, is independent of the size of the system. This differs from what it is obtained in regular lattices simulating non-mean field fracture models: there a bigger system implies a lower critical stress .

During the completion of this work, we have become aware of a similar study by Watts . He presents an interesting analytical approach to this kind of problems and solves the small-world case. Both models and the results reported, however, are different. In our model the initial distribution of intact and broken nodes after the first casualty have no contrains as in Watts’ model. Besides, and more important, the relationship between heterogeneity level, the stability of the system and the precursory activity is not accounted for in . We refer the reader to  for more details.

To sum up, we have introduced a model that accounts for the cascading events observed in many complex networks. By imposing an external pressure over the system, several magnitudes have been recorded and the system has been shown to exhibit a sort of critical point. The results point out that in order to prevent the breakdown of scale-free networks, one has to find an optimal criterion that takes into account two factors: the robustness of the system itself under repeated failures and the possibility of knowing in advance that the collapse of the system is approaching.

Y. M. thanks A. Vespignani for many useful and stimulating discussions. This work was partially supported by the Spanish DGICYT Project PB98-1594.


text:
# First weak lensing results from the Red-Sequence Cluster Survey

## 1 Introduction

In these proceedings we present some of the first weak lensing results based on the $`R_C`$band imaging data from the Red-Sequence Cluster Survey (RCS)aaahttp://www.astro.utoronto.ca/$`\stackrel{~}{}`$gladders/RCS (e.g., Gladders & Yee 2000). The Red-Sequence Cluster Survey is the largest area, moderately deep imaging survey ever undertaken on 4m class telescopes. The survey comprises 100 square degrees of imaging in 2 filters (22 widely separated patches imaged in $`R_C`$ and $`z^{}`$). Ten patches have been observed using the CFHT 12k mosaic camera, and the remaining 12 southern patches have been observed using CTIO. The depth of the survey (2 magnitudes past $`M^{}`$ at $`z=1`$) is sufficient to find a large number of galaxy clusters to $`z1.4`$.

The survey allows a variety of studies, such as constraining cosmological parameters from the measurement of the evolution of the number density of galaxy clusters as a function of mass and redshift, and studies of the evolution of cluster galaxies, blue fraction, etc. at redshifts for which very limited data are available at present.

The data are also useful for a range of lensing studies. Strong lensing by clusters of galaxies allows a detailed study of their core mass distribution. In addition, given the shallowness of the survey, the arcs are sufficiently bright to be followed up spectroscopically (e.g., Gladders, Yee, & Ellingson 2001). Thanks to the magnification of the arcs, this allows unprecedented studies of the properties of high redshift galaxies. Furthermore, in combination with detailed modeling of the cluster mass distribution, the geometry of the images can be used to constrain $`\mathrm{\Omega }_m`$.

Here we concentrate on some of the weak lensing applications, for which we use the $`R_C`$band survey data. For the purpose of weak lensing, we have to date analysed 16.4 deg2 of CFHT and 7.6 deg2 of CTIO data. The data will be used to study the mass distribution of the galaxy clusters discovered in the survey. A first study of galaxy biasing as a function of scale is presented in Hoekstra, Yee, & Gladders (2001a). In these proceedings we present some early results of the measurement of the lensing signal by large scale structure (cosmic shear) and a study of the dark matter properties of (field) galaxies. The data analysed so far do not cover complete patches, and therefore we limit the analysis to the individual pointings.

A detailed description of the data and weak lensing analysis will be provided in Hoekstra et al. (2001b; 2001c). The results presented in Hoekstra et al. (2001b) indicate that the object analysis, and the necessary corrections for observational distortions work well, which allows us to obtain accurate measurements of the weak lensing signal.

## 2 Measurement of Cosmic Shear

The weak distortions of the images of distant galaxies by intervening matter provide an important tool to study the projected mass distribution in the universe and constrain cosmological parameters (e.g., van Waerbeke et al. 2001).

The advantage of our survey, compared to other studies, is that our data have been acquired using two different telescopes, but have been analysed uniformly. The results from the two telescopes are in good agreement, suggesting that the corrections for the various observational distortions have worked well.

The left panel in Figure 1 shows the measurement of the top-hat smoothed variance for the RCS data analysed so far (16.4 deg2 from CFHT and 7.6 deg2 from CTIO). The signal-to-noise ratio of our measurements is very good, reaching $`6`$ at a radius of 2.5 arcminutes.

We use the photometric redshift distribution inferred from the Hubble Deep Field North and South to compare the observed lensing signal to CDM predictions. This redshift distribution works well, as was demonstrated by Hoekstra, Franx & Kuijken (2000) for the $`z=0.83`$ cluster of galaxies MS 1054-03. The right panel in Figure 1 shows the inferred likelihood contours for a $`\mathrm{\Lambda }`$CDM cosmology, with $`h=0.7`$. For an $`\mathrm{\Omega }_m=0.3`$ flat model we obtain $`\sigma _8=0.79\pm 0.08`$, in excellent agreement with the measurements of van Waerbeke et al. (2001).

## 3 Galaxy-galaxy lensing

Weak lensing is also an important tool to study the dark matter halos of field (spiral) galaxies (e.g., Brainerd, Blandford, & Smail 1996; Fischer et al. 2000). Rotation curves of spiral galaxies have provided important evidence for the existence of dark matter halos, but are confined to the inner regions, as are strong lensing studies. The weak lensing signal, however, can be measured out to large projected distances, and it provides a powerful probe of the gravitational potential at large radii. In particular constraints can be placed on the extent of the dark matter halos.

The lensing signal induced by an individual galaxy is too low to be detected, and one has to study the ensemble averaged signal around a large number of lenses.

For the analysis presented here, we use the 16.4 deg2 of analysed CFHT data, and select a sample of lenses and sources on the basis of their apparent $`R_C`$ magnitude. We use galaxies with $`19.5<R_C<21`$ as lenses, and galaxies with $`21.5<R_C<24`$ as sources which are used to measure the lensing signal. This selection yields a sample of 36226 lenses and $`6\times 10^5`$ sources.

The redshift distribution of the lenses is known spectroscopically from the CNOC2 field galaxy redshift survey (e.g., Yee et al. 2000), and for the source redshift distribution we again use the photmetric redshift distribution from the HDF North and South. The adopted redshift distributions give a median redshift $`z=0.35`$ for the lens galaxies, and $`z=0.53`$ for the source galaxies.

The ensemble averaged tangential distortion around galaxies with $`19.5<R_C<21`$ is presented in the left panel of Figure 2. The solid line corresponds to the best fit SIS model, for which we find an Einstein radius $`r_E=0.184\pm 0.011`$ arsec.

A better way to study the lensing signal is to compare the predicted shear (both components) from a parameterised mass model to the data. We use the truncated halo model from Schneider & Rix (1997). The results are presented in the right panel of Figure 2. With the adopted redshift distribution we obtain $`\sigma ^2^{1/2}=111\pm 5`$ km/s. It turns out that the quoted value is close to that of an $`L_{}`$ galaxy, and our results are in fair agreement with other estimates.

In addition, for the first time, the average extent of the dark matter halo has been measured. Under the assumption that all halos have the same truncation parameter, we find a 99.7% confidence upper limit of $`s<470h^1`$ kpc. More realistic scaling relations for $`s`$ give lower values for the physical scale of $`s`$ (Hoekstra, Yee, & Gladders 2001c), and therefore the result presented here can be interpreted as a robust upper limit.

We note that the results indicate a steepening of the tangential shear profile at large radii (which we interpret as a truncation of the halo), and alternative theories of gravity, such as MOND need to reproduce this. A plausible description of lensing by MOND was put forward by Mortlock & Turner (2001), which suggests that the shear in MOND drops $`1/r`$. A preliminary analysis excludes this prediction at high confidence.

## References