text:
# Mixing and condensation in a wet granular medium

## Abstract

We have studied the effect of small amounts of added liquid on the dynamic behavior of a granular system consisting of a mixture of glass beads of two different sizes. Segregation of the large beads to the top of the sample is found to depend in a nontrivial way on the liquid content. A transition to viscoplastic behavior occurs at a critical liquid content, which depends upon the bead size. We show that this transition can be interpreted as a condensation due to the hysteretic liquid bridge forces connecting the beads, and provide the corresponding phase diagram.

As is generally known from everyday experience, the mechanical properties of a granular material change dramatically if some liquid is added. The main reason is the internal cohesion due to capillary forces arising from liquid bridges between the grains Mikami et al. (1998); Israelachvili (1991); Hornbaker et al. (1997); Halsey and Levine (1998); Bocquet et al. (1998); Fraysse et al. (1999). While recent years have seen considerable progress in understanding the dynamics of dry granular materials Ristow (2000); Jaeger et al. (1996); Kadanoff (1999), the physical mechanisms underlying the properties of wet systems remain largely obscure. This is in part due to the inherent complexity of these highly nonlinear systems, but also to the fact that only few experiments have yet been conducted on this topic.

Some experiments have been recently done on the influence of liquid on the mixing behavior of systems consisting of large and small spheres Samadani and Kudrolli (2000, 2001). It was observed that as liquid was added, the tendency to segregate, which is well known for dry granulates Ottino and Khakhar (2000); Trujillo and Herrmann (February 26, 2002), was diminished. At a certain wetness, complete mixing occurred. In these experiments the morphology of heaps formed by pouring mixtures of particles with two different sizes through an orifice was investigated. In order to have a closer look at the observed effects, we have carried out experiments in which the agitation of the material is provided in a more homogeneous way, thus ensuring a more or less well defined granular temperature. Furthermore, the propensity of the system to segregate into large and small beads can be determined quantitatively in our setup. We observe a strong dependence of mixing on the liquid content, and the transition to a viscoplastic regime due to condensation of the grains into clusters, which was reported in ref.Samadani and Kudrolli (2001), was clearly observed. However, the overall behavior was found to be more complicated than previously reported, and we are able to explain all of the observed features within a rather simple theoretical model.

A typical sample consisted of $`150ml`$ of glass beads with $`R=2.5mm`$ radius, sealed in a cylindrical jar (diameter: $`80`$ mm, cf. fig. 1) together with the same volume of smaller glass beads with radius $`r`$, which was varied from $`50`$ to $`500\mu m`$, and a well defined amount of liquid. In order to prevent crystallization, the small beads were chosen to be slightly polydisperse: the spread in bead size ranged from 10 to 20 $`\%`$. The jars were mounted on a shaker, with the cylinder axis upright. The shaking motion consisted of small horizontal circles with a diameter of $`5`$ mm and a frequency of $`20`$ revolutions per second, which corresponds to a velocity of the container walls of 0.38 m/s. When the segregation pattern did not change any more (typically after one half hour), photographs of the jars were taken and analyzed. Two examples are shown in Figure 1.

We defined the mixture zone as the zone in which both the smaller and the bigger particles were visibly present. This zone was observed to be rather homogeneous concerning the mixing ratio. Occasional analysis of horizontal slices of the samples revealed a small radial variation of the mixing ratio, which is not captured by the photographs, but the results given below are qualitatively valid for the whole volume. For the results presented here we used water as the liquid, but similar results were obtained using ethanol. We can therefore exclude that leaching of ions from the glass, and other side effects due to the specific properties of water, play a significant role Olivi-Tran et al. (2002). The quantity we will discuss is the height of the mixture zone, divided by the total height of the sample, and is denoted by $`h`$. For perfect segregation, $`h0`$, while $`h=1`$ for perfect mixing. The amount of liquid added to the samples is denoted by the dimensionless quantity $`W`$, which is defined as the volume of the liquid divided by the total volume of the small beads (including the space between the beads).

In contrast to what one might anticipate on the basis of earlier work Samadani and Kudrolli (2000, 2001), the propensity of mixing of the particles did not improve monotonically with the amount of liquid added. Instead, we found in most cases that the height of the mixture zone, $`h`$, was first reduced upon addition of water. Only after a certain amount of liquid had been added did the mixing improve. In Fig. 2, $`h`$ is plotted as a function of $`W`$ for two different $`r`$. The size of the large beads was kept constant. Clearly, there is a rather well defined liquid content above which the height of the mixing zone increases when more liquid is added.

The boundary between these two regimes (i.e., the minima of the curves as in fig. 2) is indicated by the open squares in Figure 3, which represents a phase diagram of the system in the plane spanned by the radius of the small beads, $`r`$, and the liquid content, $`W`$. The dotted line serves as a guide to the eye. The full circles correspond to the occurrence of complete mixing (i.e., $`h=1`$). We can thus clearly distinguish three regimes, which we denote as gaseous, intermediate, and viscoplastic Tegzes et al. (1999, 2002), respectively. In what follows, we will discuss these regimes in some detail.

Let us start with the transition to the viscoplastic regime, in which perfect mixing is observed. This transition, which occurs as $`h1`$, appears quite abrupt when inspecting the experiment by eye. While the general motion of the whole sample was a kind of circular convection around the axis of the jar in the gaseous and intermediate regimes, the motion in the case of perfect mixing was reminiscent of kneading raisin dough. We identify this transition with the one reported recently by other authors Samadani and Kudrolli (2000, 2001). We assume that this transition into the viscoplastic regime occurs when the granular temperature of the system, i.e., the average kinetic energy per bead, is not anymore sufficient to break the liquid bridges between the beads, and is thus a kind of condensation within the granular gas, mediated by the liquid bridges.

In order to support this view, we have to discuss the involved forces in some detail. The force exerted upon two adjacent beads by the liquid bridge can be written as $`F(d,w)=2\pi r\gamma f(w)g(d)`$, where $`\gamma `$ is the surface tension of the liquid, $`r`$ is the radius of the spheres, $`w`$ is the liquid volume of the bridge, and $`f`$ is some function depending on the roughness of the sphere surfaces Halsey and Levine (1998); Willett et al. (2000). $`g(d)`$ describes the dependence of the force on the mutual distance, $`d`$, of two ideally spherical beads. The energy required to destroy the bridge is given by

$$E_B(w)=2\pi r\gamma f(w)_0^{d_p}g(d)𝑑d$$
(1)
where $`d_p`$ is the distance at which the bridge pinches off. The integral can be well approximated as $`w^{1/3}/3`$ Willett et al. (2000).

For perfect spheres and complete wetting, $`f(w)=1w`$, but for a rough surface, it tends to zero as $`w0`$. The shape of $`f(w)`$ has been discussed in detail before Halsey and Levine (1998). Based on these considerations, it can be well approximated by

$$f(w)=\left[\frac{w}{w+w_0}\right]^\mu $$
(2)
with $`w_0:=r\delta ^2`$ characterizing the roughness amplitude, $`\delta `$. For the exponent, we have $`\mu =(2\chi )/(2+\chi )`$, where $`\chi `$ is the roughness exponent of the bead surface. For (typical) scratch and dig roughness, $`\chi `$ is close to zero, and $`\mu 1`$. This corresponds to the regime found by Hornbaker et al. Hornbaker et al. (1997). In the opposite case, when $`\chi 1`$, we have $`\mu 1/3`$. We see that $`\mu `$ is confined to a narrow interval, and is very probably close to one in our system.

We now must calculate the average bridge volume, $`w`$, from the dimensionless liquid volume, $`W`$. At equilibrium, assuming that all of the liquid is in the bridges, $`w=\frac{8\pi r^3W}{3c\rho _p}`$, where $`c`$ is the coordination number of the network of bridges ($`c6`$) and $`\rho _p`$ is the packing density of the spheres ($`\rho _p0.65`$). Under dynamic conditions the typical bridge volume would be expected to be somewhat smaller than the equilibrium value. Our measurements 111M. M. Kohonen, D. Geromichalos, and S. Herminghaus; manuscript in preparation suggest that $`w=\alpha r^3W`$, where $`\alpha 0.25`$. If we now equate the bridge energy, $`E_B`$, to the average kinetic energy of the small beads, $`E_K:=\frac{m}{2}v^2`$, we obtain the implicit equation

$$r^2l^{1/\mu }\alpha ^{1/3\mu }W^{1/3\mu }r^{(21/\mu )}+\frac{\delta ^2}{\alpha W}=0$$
(3)
where we have introduced the characteristic length scale $`l:=\frac{\gamma }{\rho v^2}`$ ($`\rho `$ is the mass density of the bead material. In our case, $`\rho =2500\mathrm{k}\mathrm{g}/\mathrm{m}^3`$). Eq. (3) can be easily solved analytically for $`\mu =1`$ and $`\mu =\frac{1}{3}`$. For the most relevant case of $`\mu =1`$, we obtain

$$r=\frac{l}{2}\alpha ^{1/3}W^{1/3}\left(1+\sqrt{1\left(\frac{W_c}{W}\right)^{5/3}}\right)$$
(4)
which intersects the W-axis at the critical moisture content $`W_c=\frac{1}{\alpha }(2\delta /l)^{6/5}`$, which only depends on the roughness amplitude, $`\delta `$, but not on $`r`$. This function represents the phase boundary at which condensation sets in, and is indicated as the full curve in Figure 3. From this best-fit curve we obtain $`\delta 500nm`$, which is similar to the peak-to-peak roughness obtained from the inspection of the beads by atomic force microscopy. From the fitted value of the intrinsic length scale $`l`$, and from the known density of the glass and the surface tension of the water, we get the average granular velocity, $`v:=\sqrt{v^2}=0.11m/s`$. This compares favorably with the velocity of the shaking motion, which is expected to be of the same order of magnitude.

Let us now turn to the gaseous regime. If there is no liquid in the sample at all, segregation is observed as usual, which has become known as the ’brazil nut effect’ (BNE) Williams (1976); Rosato et al. (1987); Knight et al. (1993); M$`\ddot{\text{o}}`$bius et al. (2001). As it has been shown recently Trujillo and Herrmann (February 26, 2002), this may be understood in a continuum picture, as follows. A granular medium is strongly dissipative due to the inelastic character of the impact kinetics of the grains, which is characterized by the restitution coefficient, $`\epsilon `$. The latter is defined as the ratio of the particle velocities after and before impact, $`\epsilon =v_{rel}^{}/v_{rel}`$. Thereby the granular temperature bath, which typically corresponds to thermal temperatures in the Terakelvin range, is internally coupled to the room-temperature bath provided by the atoms constituting the grains. This inherent state of non-equilibrium is the basic reason for the peculiar features of granular matter dynamics in general.

Now consider a spherical region with radius $`R`$ within the granular gas. When a grain, or small bead, enters this region, it will have many impacts, and thereby lose a lot of energy, before it can leave this region again. When the latter is replaced by a solid sphere (i.e., a large bead of radius $`R`$), the grain will experience only one such impact. Consequently, the granular gas is hotter close to an intruder than far away in the bulk, because the intruder is a less effective heat sink. The granular gas is therefore also less dense around the intruder, and the latter rises due to the buoyancy of the ’bubble’ of reduced density it carries around itself Trujillo and Herrmann (February 26, 2002).

The height of the mixing zone is then a direct consequence of the balance between the hydrostatic pressure exerted by the large-beads fraction upon its lower interface, and the buoyancy pressure exerted upon the same interface from below by the mixed phase. This picture is supported by the observation that the large beads acquire a rather uniform spacing within the mixed phase, which might then correspond to (twice) the thickness of the reduced density layer around each large bead Trujillo and Herrmann (February 26, 2002).

The BNE increases as the energy loss upon impact is increased. This, however, is just what is achieved when liquid is added. The interparticle force exerted by the liquid has a hysteretic character Willett et al. (2000); Crassous et al. (1994); Zitzler et al. (2002), 222M. Schulz, B. M. Schulz, S. Herminghaus, Shear-induced solid-fluid transition in a wet granular medium, submitted to Phys. Rev. E and is thus intrinsically dissipative. We thus expect the restitution coefficient to be reduced, and the BNE to become more prominent, when liquid is added. More precisely, the buoyancy of the reduced density bubbles around the large beads increases, and a shallower mixing zone is needed to balance the hydrostatic pressure of the large-beads column above it. Thus the height of the mixing zone is expected to decrease when liquid is added, which is what we observe. We thus assume that throughout this regime, which is defined by a monotonically decreasing height of the mixing zone, the physics of two-particle impact is affected by the liquid merely in the sense of a decrease in the restitution coefficient, $`\epsilon `$.

Although a quantitative prediction of the BNE is yet difficult, we can give a rough estimate of the expected size of the effect. The inelastic character of the impact of the dry particles reflects itself in the BNE (and thus the tendency of demixing) of the dry material. Consequently, we expect the contribution due to the liquid to be of the same order of magnitude when the energy stored in the liquid bridges is similar to the energy lost upon impact, i.e., when $`E_B(1\epsilon ^2)E_K`$. Since $`lE_K^1`$, the corresponding moisture content, $`W_{BNE}`$, is related to $`W_c`$ via $`W_{BNE}(1\epsilon ^2)^{6/5}W_c`$. Exploiting that $`W_{BNE}`$ roughly corresponds to the open symbols in Fig. 3, we can estimate the restitution coefficient of our system. We obtain $`\epsilon =0.95\pm 0.03`$, which compares favorably with values known from the literature Ristow (2000).

Intermediate regime: As more liquid is added, the width of the hysteresis loop becomes larger, such that the kinetic energy of the impact may sometimes not be sufficient to break a bridge. It thus appears likely that small clusters will form which have a certain lifetime. Since these are heavier than the single beads, their velocity after an impact is less affected by the liquid bridges. Furthermore, the interior of a (larger) cluster may act as a sink for the liquid due to capillarity, resulting in a somewhat reduced wetness of the beads outside the cluster. Hence the ’effective’ restitution coefficient, averaged over all entities in the system, increases again. As a consequence, the BNE is now reduced, and the phases become better miscible, as we observe.

This model is in accordance with the observation that the overall appearance of the material is not expected to change noticeably in this regime: as more liquid is added, the average cluster size is just increased. The clusters take up the excess liquid, due to the increased average coordination number of the bridge network within a cluster. In other words, excess liquid is stored in liquid bridges within clusters, which are not active in the dissipation process.

We thus have elucidated the interplay between segregation (BNE) and a phase transition to viscoplasticity, and have given a robust criterion for the latter. The phase boundary can be understood from simple energy considerations.

Inspiring discussions with H. Herrmann, S. Luding, and S. Linz are gratefully acknowledged. The authors thank the German Science Foundation for financial support within the Priority Program ‘Wetting and Structure Formation at Interfaces’. MMK acknowledges the Alexander von Humboldt Foundation for generous funding.


format: text

text:
# Free–Free Absorption and the Unified Scheme

## 1 Introduction

It is an important and controversial issue what causes the low-frequency cutoff in the radio spectrum of GPS sources; synchrotron self-absorption (SSA) or free–free absorption (FFA). The discovery of FFA towards the GPS galaxy OQ 208 (Kameno et al. 2000) propounded a question how general is FFA towards GPS sources. The cold dense FFA plasma around the lobes of GPS sources could be a cocoon which smothers expansion of jets and lobes (Bicknell et al. 1997).

Such an external absorber can be a probe of the viewing angle of jets, so that a test for the unified scheme between quasars and radio galaxies (Barthel 1989), or the unified model between Seyfert 1 and 2 galaxies (Antonucci and Miller 1985) can be done. Both models presume that these classes are intrinsically identical, and suggest that the apparent differences are due to the viewing angle. With respect to Barthel’s unified scheme, he showed that the projected distances of double lobes in quasars are significantly smaller than those of radio galaxies, and proposed that a smaller viewing angle causes higher luminosities and apparent presence of the broad line components. On the Seyfert unification model, optical polarimetric observations for the Seyfert 2 galaxy NGC 1068 detected highly polarized continuum and broad Balmer lines (Miller and Antonucci 1983). Antonucci and Miller (1985) suggested that hidden broad line components arose via scattering. Presence of the broad line component, which is directly seen in Seyfert 1 but unseen in Seyfert 2 galaxies, was understood as intrinsic identity of the Seyfert 2 with Seyfert 1 galaxies. Like the unification between quasars and radio galaxies, the viewing angle was thought to play a role in different appearance.

When the jet axis is close to the line of sight, the path length through the external absorber will be longer towards the receding jet than towards the approaching jet. Thus, the FFA opacity towards a double lobe is expected to be rather asymmetric; FFA opacity should be deeper towards the receding jet. When the jet axis is nearly perpendicular to the line of sight, on the contrary, we expect rather symmetric FFA opacity towards a double lobe.

Based on this idea, we conducted a trichromatic GPS survey at 2.3, 8.4, and 15.4 GHz using the VLBA for 9 objects.

## 2 The Sample

We selected 9 sample objects, based on the GPS catalog by de Vries et al. (1997), under the criteria:

1. The peak frequency $`\nu _\mathrm{m}`$ should stand within our observing range, i.e., $`1.6`$ GHz $`<\nu _\mathrm{m}<15`$ GHz. This condition is necessary to discriminate between SSA and FFA by spectral fitting, and then to obtain the distribution of opacities.
2. Since all sources must be bright enough to be detected with the VSOP and the VLBA, we put the criteria $`S_{1.6\mathrm{GHz}}>0.1`$ Jy, $`S_{5\mathrm{G}\mathrm{H}\mathrm{z}}>0.5`$ Jy, and $`S_{15\mathrm{G}\mathrm{H}\mathrm{z}}>0.2`$ Jy.

The selected GPS sources are listed in table 1. The optical identifications are based on the NASA/IPAC Extragalactic Database (NED). Hereafter, quasars and Seyfert 1 galaxies are categorized as type-1 sources, while radio galaxies and Seyfert 2 galaxies belong to type-2 sources.

## 3 Observations and Results

The VLBA observations have been carried out on December 15, 1998. Table 1 lists the sample objects of the VLBA observations. Every object was observed at 3 frequencies with 2 to 11 scans, where each scan corresponds to integration of 11 minutes. Observations at three frequencies have been carried out almost simultaneously. The subreflector switched between the dual-frequency 2.3/8.4 GHz reflector system and the 15.4 GHz feed horn within a 1-minute gap between the cycles. We used 4 channels of 8-MHz bandwidth at 15.4 GHz, and allocated 2 channels at both 2.3 and 8.4 GHz. The correlation process was accomplished by the VLBA correlator. We applied fringe-fitting, data flagging, and a priori amplitude calibration in the NRAO AIPS. Imaging and self-calibration processes were carried out by Difmap. Synthesized beam sizes and image qualities are listed in table 1.

While relative gain errors among the antennae are corrected though amplitude self-calibration processes, further flux calibration is necessary to obtain certain spectra across the observing frequencies. For the purpose of absolute flux calibration, we also imaged four calibrators; the BL Lacertae, DA 193, 3C 279, and OT 081. These calibrators are so compact that the total CLEANed flux densities should be the same with the total flux densities measured by a reliable single-dish or a short-baseline-interferometer. Basing on the comparison between the total flux measurements by the University of Michigan Radio Astronomy Observatory (UMRAO) and the NRAO Green Bank interferometer, and summation of CLEANed flux densities, we applied flux scaling for the final results. At 15.4 GHz, for instance, total CLEANed flux densities before absolute correction were 3.276, 5.012, 25.769, and 4.172 Jy for BL Lac, DA 193, 3C 279, and OT 081, respectively. The UMRAO database provides the total flux densities at 14.5 GHz, $`3.458\pm 0.015`$, $`4.763\pm 0.035`$, $`27.69\pm 0.18`$, and $`4.380\pm 0.030`$ Jy for each, as averaged over 2 months centered at our observation. Then we calculated a correction factor ($`C=\frac{S_{\mathrm{cc}}}{S_\mathrm{U}}`$), where $`S_{\mathrm{cc}}`$ is the total CLEANed flux density of the calibrator images, and $`S_\mathrm{U}`$ is the mean flux densities of calibrators measured by the UMRAO. We have $`C=0.959\pm 0.018`$ at 15 GHz, so that we scaled our image by the factor and obtained the accuracy of the flux scale of 1.8 %. The accuracy was derived from a standard deviation of the correction factors by the four calibrators. In the same manner, we scaled images at 2.3 and 8.4 GHz by the correction factors of 0.84 and 0.79, respectively. The estimated amplitude accuracies are 4.0, 7.2 and 1.8 % at 2.3, 8.4, and 15.4 GHz, respectively.

All images with uniform weighting are shown in figure 1. Based on the images, we estimated a flux density of each component using three different methods; tasks ‘IMFIT’ and ‘JMFIT’ to adopting elliptical Gaussian distribution of the brightness, and ‘TVSTAT’ to integrate brightness within a specified area. When the resolution is insufficient to isolate a double structure, in the cases of 2.3 GHz images of 0738+313, 1333+459, and 2149+056, we applied double elliptical Gaussian fits for ‘IMFIT’ and ’JMFIT’. In these cased, we also tried Gaussian model fits in visibilities using ‘modelfit’ in Difmap instead of ‘TVSTAT’, because these images are simple enough to be fitted in visibilities while it is difficult to set adequate integral area in these images for ‘TVSTAT’. The flux densities of components are shown in table 2. Errors in the flux densities are estimated by the root sum squared of the amplitude calibration errors and the uncertainties in measurement of flux densities from the images. The uncertainties in measurement are evaluated by the standard deviation between results of the three methods, i.e. ‘IMFIT’, ‘JMFIT’, and ‘TVSTAT’ (or ‘modelfit’). Components which are considered to be lobes are labeled as A and B, in the sense that a component with a larger FFA opacity is labeled as A. The identification of lobes or a core is discussed in the next section.

Since we have observed at only three frequencies, it is impossible to discriminate between SSA and FFA by the spectral fit alone. Instead, we started from the simple FFA model assuming that SSA is negligible, and verify the assumed condition later. It should be noticed that a solid model fit to discriminate between SSA and FFA requires at least five frequencies. We are conducting VSOP observations at 1.6 and 4.8 GHz for the supplement, which will be reported in continuation. Nevertheless, global properties among type-1 and -2 GPS sources can be discussed as shown later.

Anyway, we attempted the model: an optically thin synchrotron emission from the lobe is absorbed by external FFA plasma,

$`S_\nu =S_0\nu ^{\alpha _0}\mathrm{exp}(\tau _\mathrm{f}\nu ^{2.1}).`$ (1)
Here, $`S_\nu `$ is the observed flux density in Jy, $`S_0`$ is the intrinsic flux density in Jy at 1 GHz, $`\nu `$ is the frequency in GHz, $`\alpha _0`$ is the intrinsic spectral index of the synchrotron emission, and $`\tau _\mathrm{f}`$ is the FFA coefficient. Free parameters are $`S_0`$ and $`\tau _\mathrm{f}`$ for each Gaussian component, and common $`\alpha _0`$ for all components. Therefore, the number of free parameters $`N_{\mathrm{param}}=2n+1`$, where $`n`$ is the number of Gaussian components. The number of data points is given by three frequencies times the number of Gaussian components, hence $`3n`$. Consequently, the degree of freedom will be $`3nN_{\mathrm{param}}=n1`$. When we have two Gaussian components, the degree of freedom is 1. Again, the number of frequencies is too few to verify a statistical confidence for the spectral fit. Despite this condition, the flux densities, the opacity coefficients, and the residuals of spectral fits are listed in table 3. Note that the FFA coefficient $`\tau _\mathrm{f}`$ directly corresponds to the spectral peak frequency $`\nu _\mathrm{m}`$ because the spectral peak appears at the frequency where the optical depth approximates to unity. Thus, we have

$`\nu _\mathrm{m}\tau _\mathrm{f}^{\frac{1}{2.1}}.`$ (2)
Derived best-fit spectra are shown in figures 2 and 3 for type-1 and -2 sources, respectively. Figure 4 shows the histograms of FFA opacities towards each component.

## 4 Discussion

A significant difference between type-1 and -2 GPS sources can be seen in the FFA fitting. Spectral peak frequencies of double components in type-1 sources tend to differ significantly, while those in type-2 sources are relatively similar. In other words, type-1 sources tend to show asymmetric FFA opacities towards double lobes, while type-2 sources have rather symmetric opacities, since the $`\nu _\mathrm{m}`$ and $`\tau _\mathrm{f}`$ are related by equation 2. To evaluate the asymmetry in the FFA opacities, we define the FFA opacity ratio $`R=\tau _{\mathrm{fA}}/\tau _{\mathrm{fB}}`$ ($`\tau _{\mathrm{fA}}>\tau _{\mathrm{fB}}`$). $`R`$ is an index of asymmetry in the peak frequencies $`\nu _\mathrm{m}`$, too, as

$`R=\left({\displaystyle \frac{\nu _{\mathrm{mA}}}{\nu _{\mathrm{mB}}}}\right)^{2.1}.`$ (3)
Even if the spectral peak is caused by SSA, beyond our assumption, $`R`$ represents asymmetry in terms of SSA. The spectrum of power-law synchrotron radiation with SSA is expressed by

$`S_\nu =S_0\nu ^{2.5}\left[1\mathrm{exp}(\tau _\mathrm{s}\nu ^{\alpha _02.5})\right],`$ (4)
where $`\tau _s`$ is the SSA opacity coefficient. The spectral peak $`\nu _\mathrm{m}`$ approximates to $`\tau _\mathrm{s}^{\frac{1}{2.5\alpha _0}}`$. Then the relation between $`R`$ and $`\tau _\mathrm{s}`$ will be

$`R=\left({\displaystyle \frac{\tau _{\mathrm{sA}}}{\tau _{\mathrm{sB}}}}\right)^{\frac{2.1}{2.5\alpha _0}},`$ (5)
though $`R`$ is derived from the spectral fits using the FFA model.

The histogram of $`R`$ clearly exhibits the difference between type-1 and -2 sources (see figure 4). The weighted means of $`R`$ for type-1 and -2 are $`4.97\pm 1.15`$ and $`1.36\pm 0.51`$, respectively. Here, we perform a statistical test, whether any significant difference arises between the FFA opacity ratios of type-1 and -2 groups. Let us put a testing hypothesis, which assumes $`R`$ is identical for these two subsets. Then, the $`T`$ value defined as

$`T={\displaystyle \frac{|\overline{X_1}\overline{X_2}|}{\sqrt{\left(\frac{1}{n_1}+\frac{1}{n_2}\right)\left(\frac{S_1+S_2}{n_1+n_22}\right)}}},`$ (6)
must follow the $`t`$ distribution with the degrees of freedom of 7. Here, $`\overline{X_1}`$ and $`\overline{X_2}`$ are mean values of a variable $`X`$ in subsets 1 and 2, respectively, $`n_1`$ and $`n_2`$ are the number of data; $`S_1`$ and $`S_2`$ are the sum of residual squared. The calculated value $`T=3.90>t(7,0.01)=3.5`$ rules out the testing hypothesis. Thus, the opacity ratio of type-1 sources is significantly larger than that of type-2 sources, with the confidence larger than 99 %. Even if we take the T-test for logarithms $`\mathrm{ln}R`$, the values are $`1.70\pm 0.36`$ and $`0.57\pm 0.52`$ for type-1 and -2 sources, respectively. Then, $`T=4.26>t(7,0.01)`$ suggests, too, that type-1 sources are more asymmetric than type-2 sources in terms of FFA opacity.

The result is consistent with Barthel’s unified scheme between RGs and QSOs (Barthel 1989). If the line of sight is close to the jet axis, as thought to be type-1, a large difference of the path length in external plasma towards the lobes causes an asymmetric FFA. In case of type-2 sources, the line of sight is nearly perpendicular to the jet axis, so that a small difference in the path length results in a relatively symmetric FFA.

To check if any bias is included in the statistics, we also test the distribution of intrinsic spectral indices $`\alpha _0`$ and ratios of intrinsic flux densities $`FR=\frac{S_{0\mathrm{A}}}{S_{0\mathrm{B}}}`$. The spectral indices are $`0.64\pm 0.16`$ and $`1.12\pm 0.25`$ for type 1 and 2, respectively. $`T=3.06`$ between $`t(7,0.05)=2.365`$ and $`t(7,0.01)=3.5`$ suggests that $`\alpha _0`$ is likely to be different between the two classes. The flux density ratios $`\mathrm{ln}FR`$ are $`1.76\pm 0.85`$ and $`0.17\pm 0.62`$ for type 1 and 2, respectively, and $`T=2.74<t(7,0.1)`$ indicate no significant difference.

The larger $`\alpha _0`$ for type-1 objects probably indicates contamination of the core component. In fact, the flat spectrum of component A in 0646+600 and 1333+459 implies that these components are the core. Alternatively, let us use $`\alpha _0`$ of component B (lobes or jets) for these two objects, i.e., $`\alpha _0=1.54`$ and $`1.45`$ for 0646+600 and 1333+459, respectively, as are derived in Case 2. The mean $`\alpha _0`$ for type-1 objects will be $`0.95\pm 0.47`$, which attributes to no significant difference between that of type-2 objects with $`T=0.59<t(7,0.1)`$.

From the above statistical considerations we conclude that the two groups show no intrinsic difference, but are apparently different in terms of FFA opacities. What does this mean? Here, we discuss the reason of the apparent difference of the opacity ratios.

Case 1: The viewing angle of the jet in type-1 source is smaller than that of type-2 sources

This idea is consistent with the unified scheme (Barthel 1989). It simply explains why the FFA opacities are asymmetric towards type 1s, while no difference in intrinsic properties between type 1 and 2 is found. The significant difference of the $`\tau _{\mathrm{fB}}`$ between type-1 and -2 objects (see figure 4) can be understood in this model: the approaching component would have less path length in the ambient FFA plasma as shown in figure 7. The external absorption model is self-consistent with the FFA fitting. This model does not require that SSA is intrinsically asymmetric towards the lobes.

Case 2: Type-1 sources consist of the core and one-sided jet

In this case, the inverted spectrum towards the core component is likely due to SSA. The intrinsic spectral index $`\alpha _0`$ of a core could be larger than that of a jet or a lobe. As shown in figure 6, spectral indices of type-1 sources are likely to be larger than those of type-2 sources. For example, the spectral fits for 0646+600 and 1333+459 will be much better if a larger $`\alpha _0`$ is put for components A. In the case of 0646+600, $`\alpha _0=0.11`$ and $`1.54`$ for components A and B result in the best fit, though we lose degrees of freedom, and we have $`S_{0\mathrm{A}}=1.19\pm 0.02`$, $`\tau _{\mathrm{fA}}=6.22\pm 0.25`$, $`S_{0\mathrm{B}}=8.32\pm 0.14`$ and $`\tau _{\mathrm{fB}}=9.91\pm 0.25`$. In the case of 1333+459, $`\alpha _0=0.50`$ and $`1.45`$ for A and B, respectively, are the best fit to derive $`S_{0\mathrm{A}}=1.01\pm 0.02`$, $`\tau _{\mathrm{fA}}=7.7\pm 0.3`$, $`S_{0\mathrm{B}}=4.89\pm 0.11`$ and $`\tau _{\mathrm{fB}}=7.8\pm 0.3`$. In both cases, components A are likely to be a core, rather than a lobe, because of their flat spectrum. If we assume intrinsic bipolarity of the jet, it is necessary to consider why the counterjet is unseen.

One may consider that the Doppler boosting effect possibly causes the apparent one-sided jet. In this interpretation, component B is approaching towards us to be amplified, while the unseen counterjet is receding to be dimmed. However, the Doppler boosting effect often results in variability of flux densities and a large polarization degree. These properties are unlikely for GPS sources.

An alternative interpretation is that the unseen counterjet is severely obscured via FFA. This could be an attractive model in which the diffuse emissions opposite to components B in 0646+600 and 0738+313 can be understood as the counterjets. Deeper imaging capabilities at higher frequencies are necessary to confirm the counterjet.

Whatever attenuates the counterjet, a smaller viewing angle in type-1 sources than that of type-2 sources is implied. Thus, this case is also consistent with the Barthel’s unified scheme. The path length in the external absorber towards component B of type-1 sources is expected to be shorter than those of type-2 sources (see figure 7). Smaller $`\tau _{\mathrm{fB}}`$ of type-1 sources than that of type-2 sources can be understood in this context.

Case 3: Type-1 sources are smaller than type-2 sources.

When the plasma density decreases as a function of the radius from the nucleus, like the isothermal King model (King 1972), denser FFA opacities near the center will be produced. The asymmetry of the opacity can be enhanced, even if the viewing angle remains the same, if the source size is small (Kameno et al. 2001). This idea is similar to the hypothesis that the FFA plasma is denser towards type-1 sources than towards type-2 sources, as is scaled by the core radius of the ambient FFA plasma. However, this model requires that the mean opacity of type-1 sources should be larger than that of type-2 sources. The histogram $`\tau _\mathrm{B}`$ (figure 4) shows the opposite behavior, so that this idea is not supported.

Consequently, the simplest explanation for the difference in the asymmetry of opacity is the difference of viewing angle. This should be discussed coupled with the absorption mechanisms, FFA or SSA. If SSA dominated in the low-frequency cutoff, we expect that the peak frequencies of components A and B would be statistically unbiased in the co-moving frame with each synchrotron emitter. When the synchrotron sources emanate from the core at a relativistic speed, the peak frequencies of approaching and receding components would become higher and lower, respectively, due to the Doppler effect. Together, the flux densities of approaching and receding components would be also amplified and dimmed, respectively, by the Doppler effect. Hence, the flux density ratio $`FR=S_{0\mathrm{A}}/S_{0\mathrm{B}}`$ is expected to correlate with the peak frequency. However, no significant correlation is found between $`R=\tau _{\mathrm{fA}}/\tau _{\mathrm{fB}}`$ which is related to the peak frequency, and $`FR`$, with the correlation coefficient of $`\rho =0.29`$ (see figure 6). These statistics of $`R`$ and $`FR`$ does not support the simple SSA-only model, and requires that double lobes are intrinsically asymmetric. Saikia et al. (1995) reported that compact steep spectrum sources with small linear size are intrinsically asymmetric. Carvalho (1998) showed that interaction of the jet with a non-homogeneous intragalactic medium can result in asymmetric evolution. Since the peak frequency in terms of SSA is anti-correlated with the linear size (O’Dea et al. 1998), the intrinsic asymmetry of the SSA opacity can be produced. The presence of SSA and its intrinsic asymmetry cannot be ruled out, nevertheless, the unified scheme coupled with the FFA model can simply account for the difference in asymmetry of opacities.

## 5 Conclusions

VLBA observations for nine sources at 2.3, 8.4, and 15.4 GHz have been carried out, to reveal the morphologies of all objects. Spectral model fitting is applied to obtain spatial distribution of FFA opacities towards the radio emission of individual sources. A difference between type-1 (quasars and Seyfert 1 galaxies) and type-2 (radio and Seyfert 2 galaxies) is found, in terms of the FFA opacity ratio $`R=\tau _{\mathrm{fA}}/\tau _{\mathrm{fB}}`$ between two components A and B of each object. Type-1 objects tend to show significantly larger opacity ratios than type-2 sources do. Asymmetry in FFA opacities suggests that path lengths through ambient absorbers towards twin lobes are significantly different. Therefore, larger opacity ratios of type 1 objects indicate the axes between lobes are smaller than those of type 2 objects (see figure 7). This result supports the unified scheme between quasars and radio galaxies proposed by Barthel (1989).

## Acknowledgments

We the VLBA operated by the National Radio Astronomy Observatory, which is a facility of the National Science Foundation (NSF). This work was supported by Grant in Aid (c) 13640248 from Japan Society for the Promotion of Science.

## References

Antonucci, R. R. J., & Miller, J. S. 1985, ApJ, 197, 621

Barthel, P. D. 1989, ApJ, 336, 606

Bicknell, G. V., Dopita, M. A., & O’Dea, C. P. 1997, ApJ, 485, 112

Carvalho, J. C. 1998, A&A, 329, 845

Kameno, S, Horiuchi, S., Shen, Z.-Q., Inoue, M., Kobayashi, H., Hirabayashi, H., & Murata, Y. 2000, PASJ, 52, 209

King, I. R. 1972, ApJ. 174, L123

Kameno, S., Sawada-Satoh, S., Inoue, M., Shen, Z.-Q., & Wajima, K. 2001, PASJ, 53, 169

de Vries, W. H., Barthel, P. D., & O’Dea, C. P. 1997, A&A, 321, 105

Miller, J. S., & Antonucci, R. R. J. 1983, ApJ, 271, L7

O’Dea C. P. 1998, PASP, 110, 493


text:
# Distance, dissimilarity index, and network community structure

## I Introduction

A graph (network) of vertices (nodes) and edges is a useful tool in describing the interactions between different agents of a complex system. For example if we want to analyze protein-protein physical interactions in yeast Saccharomyces cerevisiae Uetz et al. (2000), we would like to denote each protein as a distinct vertex of a graph, and setup an edge between two vertices if the corresponding proteins have direct physical interactions. Many such kinds of networks are constructed in sociological, biological, and technological fields, and they usually have very complicated connection patterns. What one needs is a method that is capable of classifying vertices of a complex network into different clusters (communities). If a network is appropriately decomposed into a series of functional units, (a) the structure of the network can be better understood and the relationship between its different components will be clear, (b) the principal function of each cluster can be inferred from the functions of its members, and (c) possible functions for members of a cluster can be suggested by comparing the functions of other members. Network clustering techniques are therefore very important in the emerging fields of bioinformatics and proteomics.

A good clustering method needs to satisfy two conditions: First, the inherent structure of the network should be reserved; Second, it should provide a quantified resolution parameter to mark the significance of the clusters obtained at each level of the partitioning process. The global organization of a network should already be identified at low resolutions and more and more fine structures emerge as the resolving power is increased.

Many existing methods Wasserman and Faust (1994); Ravasz et al. (2002) only take account of local information of each vertex, such as number of nearest-neighbors shared with other vertices, number of vertex-independent paths to other vertices, etc.. Recently, Girvan and Newman Girvan and Newman (2002) suggested an elegant global algorithm which extended the concept of vertex betweenness centrality of Freeman Freeman (1977) also to edges. Their algorithm works iteratively by removing the current edge(s) of the highest degree of betweenness centrality. When applying to an ensemble of random modular networks, this algorithm greatly outperforms some conventional methods Girvan and Newman (2002). On the other hand, it does not provide a parameter to quantify the differences between communities.

In reference zho (a) a Brownian particle is $`\mathrm{`}\mathrm{`}`$introduced” into a network to $`\mathrm{`}\mathrm{`}`$measure” the distances between vertices. In the present work, we extend the basic idea of zho (a) by defining, based on this distance matrix, a quantity called the dissimilarity index between nearest-neighboring vertices. The dissimilarity index signifies to what extent two nearest-neighboring vertices would like to be in the same community. A hierarchical algorithm is then worked out; it takes use of information on the dissimilarity indices and decompose a network into a hierarchical sequence of clusters. Each of the communities is characterized by an upper and a lower dissimilarity threshold.

The method, which could work on unweighted as well as weighted networks, is applied to several artificial and real networks, and very satisfying results are obtained. For the case of random modular networks, the present algorithm outperforms the method of Girvan and Newman Girvan and Newman (2002). When applying the algorithm to the protein-protein interaction network of yeast, we are able to identify many protein clusters which have well defined biological functions.

In section II we review the distance measure of reference zho (a) and define a dissimilarity index for each pair of nearest-neighboring vertices. A dissimilarity-index-based hierarchical algorithm is outlined in section III, and applied to two kinds of artificially generated networks and four real-world networks in section IV. We conclude our work in section V.

## II Distance measure and dissimilarity index

In the opinion of Flake, Lawrence, and Giles Flake et al. (2000), a community in a (sub)graph should satisfy the requirement that each vertex’s total intra-community interaction be stronger than the total interaction with other vertices in the (sub)graph. This turns out to be a very strong constraint. In this work, we weaken this condition and require only that a vertex should have stronger total interaction with other vertices of its own community than with vertices of any another community of the (sub)graph.

We consider a connected network of $`N`$ vertices and $`M`$ edges. The network’s connection pattern is specified by the generalized adjacency matrix $`A`$. We assume that the value of each non-zero element of matrix $`A`$ (say $`A_{ij}`$) denotes the interaction strength between vertex $`i`$ and $`j`$. The distance, $`d_{ij}`$, from vertex $`i`$ to vertex $`j`$ is defined as the average number of steps needed for a Brownian particle on this network to move from vertex $`i`$ to vertex $`j`$ zho (a). At each vertex (say $`k`$) the Brownian particle will jump in the next step to a nearest-neighboring vertex (say $`l`$) with probability $`P_{kl}=A_{kl}/_{m=1}^NA_{km}`$. The distance matrix thus defined is asymmetric (in general $`d_{ij}d_{ji}`$), and it is calculated by solving $`N`$ linear-algebraic equations zho (a).

Taking any vertex $`i`$ as the origin of the network, then the set $`\{d_{i1},\mathrm{},d_{i,i1},d_{i,i+1},\mathrm{},d_{iN}\}`$ measures how far all the other vertices are located from the origin. Therefore it is actually a perspective of the whole network with vertex $`i`$ being the viewpoint. Suppose vertex $`i`$ and $`j`$ are nearest-neighbors ($`A_{ij}>0`$), the difference in their perspectives about the network can be quantitatively measured. We define the dissimilarity index, $`\mathrm{\Lambda }(i,j)`$, by the following expression:

$$\mathrm{\Lambda }(i,j)=\frac{\sqrt{\underset{ki,j}{\overset{N}{}}[d_{ik}d_{jk}]^2}}{(N2)}.$$
(1)
If two nearest-neighboring vertices $`i`$ and $`j`$ belong to the same community, then the average distance $`d_{ik}`$ from $`i`$ to any another vertex $`k`$ ($`ki,j`$) will be quite similar to the average distance $`d_{jk}`$ from $`j`$ to $`k`$, therefore the network’s two perspectives (based on $`i`$ and $`j`$, respectively) will be quite similar. Consequently, $`\mathrm{\Lambda }(i,j)`$ will be small if $`i`$ and $`j`$ belong to the same community and large if they belong to different communities.

## III the algorithm

We exploit the dissimilarity index to decipher the community structure of a network. After the distance matrix $`\{d_{ij}\}`$ and the dissimilarity indices for all the nearest-neighboring vertices $`\{\mathrm{\Lambda }(i,j)\}`$ are obtained, the algorithm works as follows:

Intially the whole network is just one single community. This community is assigned an upper dissimilarity threshold $`\theta _{\mathrm{upp}}`$ equalling to the maximum value of all the different dissimilarity indices.

For each community, a resolution threshold parameter $`\theta `$ is introduced and is assigned the initial value $`\theta _{\mathrm{upp}}`$ of that community. The algorithm is unable to discriminate between two nearest-neighboring vertices $`i`$ and $`j`$ when $`\mathrm{\Lambda }(i,j)\theta `$; if this happens, vertices $`i`$ and $`j`$ are marked as $`\mathrm{`}\mathrm{`}`$friends”.

The $`\theta `$ value is decreased differentially. All edges in the community are examined to see whether two nearest-neighboring vertices are friends. Different friends sets of the community are then formed, each of which contains all the friends of the vertices in the set. There may also be vertices in the community that do not have any friends. Each of these vertices is moved to the friends set that has the strongest interaction with it. After this operation, vertices of the community are distributed into a number of disjointed sets (this number may be unity).

Each vertex in a subcluster should have stronger interaction with vertices within this subcluster than with vertices of any another subcluster of this community. To fulfill this requirement, we perform a local adjustment process: move each of the vertices that fail to meet this requirement to the friends set that has the strongest total interaction with it. This adjustment process is performed simultaneously for all these unstable vertices and is repeated until no unstable vertices remains.

If vertices of the community remain together, the algorithm returns to step 3. If they are divided into two or more sets, then the community under processing is assigned a lower dissimilarity threshold $`\theta _{\mathrm{low}}`$ equalling to the current $`\theta `$ value, and it is no longer considered. Each of the identified subsets of this community is regarded as a new (lower-level) community, with upper dissimilarity threshold $`\theta _{\mathrm{upp}}`$ equalling to the current $`\theta `$ value. The algorithm returns to step 2 to work with another identified community.

After all the (sub)communities are processed, a dendrogram is drawn to demonstrate the relationship between different communities as well as the upper and lower dissimilarity thresholds of each community. The vertex set of each community is also reported.

The above procedure could be easily implemented with C++ programming language. The source code as well as the data for the examples studied in the following section will be made publicly available zho (b).

## IV Applications

We test the performance of the above-mentioned algorithm by applying it first to two kinds of artificial networks and to three real-world networks.

### IV.1 Artificial random modular networks

To quantitatively compare with the work of Girvan and Newman Girvan and Newman (2002) the algorithm is first applied to a random modular network. The network has $`128`$ nodes, which are divided into $`4`$ modules of size $`32`$ each. Each vertex has on average $`16`$ edges connecting to other vertices, and on average $`\overline{z}_{\mathrm{out}}`$ of each vertex’s edges are to vertices of other modules. All the edges are setup randomly with these two fixed expectation values. The present method is able to recover the modular structure of the network up to $`\overline{z}_{\mathrm{out}}7`$. It slight outperforms the method of Girvan and Newman Girvan and Newman (2002) in performance. For example, working on an ensemble of random graphs with $`\overline{z}_{\mathrm{out}}=6.0`$ by the present method, on average only $`4.5`$ vertices are misclassified, each of which is assigned a cluster identity different from those of the majority of vertices of its module; while on average about $`13`$ vertices are misclassified by the method of Girvan and Newman Girvan and Newman (2002).

In figure 1, the community structure of a randomly generated modular network with $`\overline{z}_{\mathrm{out}}=6.0`$ is demonstrated. When the resolution threshold is beyond $`0.323`$, the network as a whole could be regarded as a giant community. At resolution threshold $`0.323`$, however, $`3`$ subgroups suddenly emerge, with size $`32`$, $`32`$, and $`64`$, respectively. The first two communities correspond to two modules of the network, and the last one is the merge of the other two modules. At resolution threshold $`0.319`$, this later community again is divided into two subcommunities of $`32`$ vertices each, corresponding to the remaining two modules. At resolution threshold $`0.258`$, one of the modules of the network is found to fission into two subgroups of size $`14`$ and $`18`$, respectively. In this example, the designed four modules of the network correspond to the resolution range from $`0.258`$ to $`0.319`$.

How to interpret the resolution parameters in the dengrograms such as that shown in figure 1? Take module $`2`$ and module $`3`$ as examples. Figure 1 suggests that edges between these two modules have dissimilarity indices larger than $`0.323`$, while edges within these modules have dissimilarity indices $`0.227`$. Therefore there is a large dissimilarity gap of about $`0.1`$ between an inter-modular edge and an intra-modular edge.

It is noticeable that by the present algorithm, each community has certain range of stability. Subcommunities emerge only when the resolution threshold is lowered below certain level, and they emerge abruptly.

### IV.2 Regular hierarchy networks

We analyze here the community structure of the model hierarchy network studied by Ravasz and coauthors Ravasz et al. (2002). The network is constructed by several steps Ravasz et al. (2002): At level $`n=0`$, a fully connected unit of four vertices is generated. At level $`n=1`$, three replicas of this unit are added and the external vertices of these replicas are connected to the central vertex of the $`n=0`$ unit, while the central vertices of the replicas are connected to each other. This replication-connection process could be continued to any desired level $`n`$. In figure 2A such a network at level $`n=2`$ is shown. It was remarked Ravasz et al. (2002) that conventional network clustering methods are unable to uncover the hierarchical structure of such a network. The present method, however, works very well: figure 2B demonstrates the obtained community structure of the network figure 2A. The hierarchy organization of the vertices in the network is largely reserved in 2B. At resolution threshold $`1.95`$ the network is divided into $`4`$ subgroups of size $`3`$ and a giant component of size $`62`$. Later at resolution threshold $`1.89`$, this giant component again is fissioned into $`2`$ parts: one part has size $`12`$ and is further divided into $`3`$ subgroups of size $`4`$ at resolution threshold $`1.52`$; the other part has size $`50`$, which, at resolution threshold $`1.53`$ further decomposes into $`3`$ subgroups of size $`13`$, $`13`$, and $`14`$, respectively. At resolution threshold $`0.91`$, each of these three subgroups is further divided into $`3`$ subgroups.

### IV.3 The karate club network

The karate club data Zachary (1977) examined in references Girvan and Newman (2002) and zho (a) is re-evaluated here. This network is weighted, each edge is assigned a different strength. The present algorithm leads to the community structure of figure 3. At resolution threshold $`1.67`$ the network decomposes into one small component of $`5`$ vertices and a large component of $`29`$ vertices. At resolution threshold $`0.87`$, this large component further decomposes into two subgroups, One of which has $`18`$ members and the other has $`11`$ members. Comparison with the actual fission pattern is also shown in figure 3.

### IV.4 The foot-ball team network

The foot-ball team network compiled by Girvan and Newman Girvan and Newman (2002) and studied in references Girvan and Newman (2002) and zho (a) is re-investigated here. The present method results in the community structure of figure 4. Each vertex’s actual group-identity is also shown for comparison. In the resolution region between $`0.41`$ and $`0.64`$ there are $`12`$ communities according to the present algorithm. Of the $`12`$ actual groups, only members from group-$`12`$ are distributed to other groups (with good reasons, because actually there are very few direct interactions between the five members of this cluster). Vertex $`111`$ are classified together with members of group-$`11`$, we have checked that this vertex has $`8`$ edges linking to group-$`11`$ and only $`3`$ edges to other groups. Vertex $`59`$ is classified together with members of group-$`9`$, we have also checked that it has stronger interaction with group-$`9`$ than with any another group.

The organization of the different teams suggested by the present algorithm seems to be even better than the original organization.

### IV.5 The scientific collaboration network

The scientific collaboration network compiled by Girvan and Newman Girvan and Newman (2002) and examined in references Girvan and Newman (2002) and zho (a) is also re-examined. This network is also weighted. The present method suggests a community structure shown in figure 5. In accordance with the actual situation, on the global scale, the network clearly has $`3`$ giant communities of comparable sizes. Each of these giant communities could further be decomposed into several subcommunities when the resolving power is increased.

### IV.6 The protein interaction network of yeast

The protein interaction network of yeast is constructed based on the data reported in references Xenarios et al. (2000) and Deane et al. (2002), it contains $`1471`$ proteins and $`2770`$ edges (protein-protein physical interactions). This network has already been studied in reference zho (a); here we constructed a reduced interaction network based on the original one. First, self-connection is removed; second, proteins which are connected to the network by only one edge are removed. The second step is continued until no proteins of degree one remains. The reason to remove all the proteins of degree one is that, according to the idea of Girvan and Newman Girvan and Newman (2002), a vertex that is connected to the network by just one edge should be in the same community as its nearest-neighboring vertex, therefore its status need not to be considered separately. Of cause, we have checked that actually identical results are obtained when the network-reduction process is not performed. The reduced network contains $`871`$ proteins and $`2043`$ unweighted interactions (edges).

The community structure of this network is demonstrated in figure 6. It seems to be strikingly different from those of the other networks studied in this paper. At the resolution range between $`1.5`$ to $`18.0`$ many small communities appear, but the network is dominated by just one large cluster of size proportional to the total size of the network. This is in accordance with reference zho (a) where the original network was decomposed into one large component and several small components. When the resolution threshold is decreased below $`1.5`$, the largest cluster is divided into several subclusters of comparable sizes. The biological significance of such a community structure is yet to be investigated.

Based on the community structure shown in figure 6, we can construct clusters of proteins that might be of biological significance. Here we just show three examples of such protein clusters, corresponding respectively to higher, medial, and lower resolution thresholds.

The first example is a cluster which appears at resolution threshold $`18.04`$. It contains $`16`$ proteins and $`33`$ edges, and has the structure shown in figure 7A. This cluster is stable, namely that each vertex in it is more connected to vertices in this cluster than to vertices outside; and it has no further subcommunity structure. According to the protein interaction databank Xenarios et al. (2000); Bairoch and Apweiler (2000), $`15`$ of these proteins are all involved in ATP synthesis process in yeast. They may form a very important part of yeast’s mitochondrial ATPase complex. One protein of this cluster, YIL124W, is a hypothetical membrane protein. Because this last protein has only one interaction with other members of the cluster, it may not have similar biological functions as the other members.

The second example is a cluster which appears at resolution threshold $`5.11`$. It contains $`11`$ proteins and $`38`$ edges, and has the structure shown in figure 7B. This cluster is also stable and has no further structure. According to the protein interaction databank Xenarios et al. (2000); Bairoch and Apweiler (2000), among these $`11`$ proteins, YBL084C, YFR036W, YHR166C, YKL022C are known to be cell division control proteins; YGL240W plays a role in cell cycle and mitosis; YDR118W, YNL172W, YOR249C probably are membrane proteins; and YLR127C, YDL008W, YLR102C are hypothetical proteins whose functions remain to be determined. It is quite likely that all the proteins in this cluster are closely involved in cell division and membrane fission process. We anticipate that the three hypothetical proteins of this cluster will also have similar biological functions.

The third example is a cluster which appears only when the resolution threshold is refined to below $`0.88`$. It contains $`14`$ proteins and $`41`$ protein-protein interactions. This cluster is also stable and has no further structure. The interaction pattern of this cluster is demonstrated in figure 7C. Among these $`14`$ proteins, according to the protein interaction databank Xenarios et al. (2000); Bairoch and Apweiler (2000), YCR093W, YPR072W, YDL165W, YER068W, YIL038C are general negative regulator of transcription subunits; YAL021C is a glucose-repressible alcohol dehydrogenase transcriptional effector; YNR052C is a ubiquitous transcription factor; YDR443C, YGR104C are suppressors of RNA polymerases; YNL025C is the RNA polymerase II holoenzyme cyclin-like subunit; YPL042C is the meiotic mRNA stability protein kinase UME5; YGR092W is the cell cycle protein kinase DBF2; and YKR036C and YFL028C are two hypothetical proteins. It is quite likely that this cluster is mainly involved in RNA transcription process and we also anticipate that the two hypothetical proteins of this cluster are strongly related with this biological function.

To conclude this subsection, we stress that, based on the community structure of figure 6 many clusters of proteins can be constructed. Here we have mentioned just three examples. These identified protein clusters could help researchers to assign possible biological functions to hypothetical proteins, and could also suggest possible proteins that may be involved in carrying out a particular biological reaction.

## V conclusion and discussion

In our earlier work zho (a), the distance between two vertices of a graph is defined as the average number of steps a Brownian particle takes to move from one vertex to the other. Based on this distance measure, in the present work we define a dissimilarity index to signify to what extent two nearest-neighboring vertices will be different from each other. We observe that vertices belonging to the same group usually have very small dissimilarity indices between them, while vertices of different communities usually have large dissimilarity indices between them. The observation leads naturally to an algorithm of network clustering. We applied this method to several artificial networks and also to different real networks in social and biological systems and satisfactory results are obtained. Different clusters of a network obtained by our method are characterized by a range of resolution threshold.

The examples studied by us in this paper suggest that our algorithm is very promising in identifying the community structure of a complex networked system. Why it works? Maybe it is because of the following reasons. First, the vertex-vertex distance measure has taken into account the topological structure of the network as well as the local connections of the network. The distances from one vertex to all the other vertices of the network actually give a perspective of the whole network viewed from this vertex. Second, the dissimilarity index defined by equation (1) compares the perspectives viewed from two nearest-neighboring vertices. It is intuitively appealing to assume that the perspectives of the different vertices of the same community are similar to each other while those of vertices of different communities will be quite different.

It is anticipated that the present work will find applications in the field of complex networks, as well as in the fields of sociological and biological sciences.

## Acknowledgement

This research is made possible by a post-doctoral fellowship from the Max-Planck Society. The author is grateful to Professor Reinhard Lipowsky for his constant support.


text:
# Thermodynamic Properties and Pressure Effect on the Superconductivity in CaAlSi and SrAlSi

## I Introduction

Since the discovery of superconductivity at 39 K in $`MgB_2`$, the attention has turned to other intermetallic compounds with similar structure or lattice symmetry. Besides $`MgB_2`$ only a few binary compounds between the IIA and IIIB elements crystallize in the hexagonal $`C32`$ structure, i.e. $`CaGa_2`$, $`SrGa_2`$, and $`BaGa_2`$. These gallium-based binary compounds are not superconducting.1  However, a new class of pseudo-ternary compounds with $`C32`$ structure, $`AX_{2x}Si_x`$ with $`A=Ca,Sr,Ba`$ and $`X=Ga,Al`$, was recently synthesized and characterized.1 ; 2 ; 3 ; 4 ; 5 ; 6  The pseudo-ternary compounds are derived from the binary ones by partially substituting $`X`$ (e.g. $`Ga`$) by silicon which results in a change of the lattice parameters1 ; 4  and, presumably, of the electronic structure. The compounds are isostructural to $`MgB_2`$ where the $`A`$-ions occupy the $`Mg`$ sites and the $`X`$/$`Si`$ ions are randomly distributed among the boron sites in the honeycomb planes of the $`C32`$ structure. The search for superconductivity was successful for most of the materials mentioned above, however, stoichiometric $`BaAlSi`$ is not superconducting above 2 K (the measuring limit of previous investigations). Magnetization measurements show the characteristic features of type II superconductors.2  The $`C32`$ structure is stable within a wide range of $`Al`$:$`Si`$ and $`Ga`$:$`Si`$ ratios and superconductivity was observed in a range from 0.6 to 2.3 ($`Al`$:$`Si`$) and 0.3 to 1.9 ($`Ga`$:$`Si`$), respectively.1 ; 4  The maximum $`T_c`$’s were found near the stoichiometric 1:1:1 composition ($`x`$=1) in $`CaAl_{2x}Si_x`$, $`SrGa_{2x}Si_x`$, and $`BaGa_{2x}Si_x`$ but at $`x`$=0.85 in $`CaGa_{2x}Si_x`$.

Because of the structural similarity to the famous sister compound, $`MgB_2`$, the question has been raised whether the properties and mechanisms leading to superconductivity in the pseudo-ternary compounds are of the same origin as those in magnesium diboride and why their $`T_c`$ is so much lower (the highest $`T_c`$ $``$ 7.9 K was found in $`CaAlSi`$). We have shown, based on thermoelectric power measurements, that the majority of charge carriers in the pseudo-ternary compounds are electron like in contrast to the $`\sigma `$-hole carriers in $`MgB_2`$.1 ; 4  This conclusion is supported by recent band structure calculations showing that the electron-like d-states of the alkaline-earth metals contribute up to 60 % to the total density of states at the Fermi energy.7

It is interesting to note that there is no obvious correlation of $`T_c`$ with the mass of the alkaline-earth ion in the compound $`A`$-$`GaSi`$. However, in $`A`$-$`AlSi`$ $`T_c`$ decreases systematically with increasing mass of the $`A`$-ion from 7.9 K ($`A=Ca`$) to 5 K ($`A=Sr`$) and $`BaAlSi`$ is not superconducting above 2 K.4 ; 6  Imai et al.6  attributed the trend of $`T_c`$ in the $`A`$-$`AlSi`$ system to a decrease of the density of states, $`N(E_F)`$, due to the increasing mass of the $`A`$-ion assuming thereby that the Debye temperature (the pre-factor in the BCS formula for $`T_c`$) does not change if $`Ca`$ is replaced by $`Sr`$ and $`Ba`$. This assumption, however, is not justified because, unlike in $`MgB_2`$, the electron-phonon interaction leading to the formation of Cooper pairs and superconductivity is expected to involve the vibrations of the alkaline-earth ions and their mass should affect the superconducting $`T_c`$. The latter argument is supported by the results of recent band structure calculations for $`A`$-$`AlSi`$ indicating that $`N(E_F)`$ actually increases if $`Ca`$ is replaced by $`Sr`$ or $`Ba`$.7  Furthermore, the alkaline-earth d-states contribute to the Fermi surface to a large extent and there is no reason to assume that they do not participate in the electron-phonon coupling. Additional investigations are needed to reveal the physical nature of the superconducting state.

We, therefore, have conducted measurements of the specific heat as well as the effect of hydrostatic pressure on the superconducting transitions of $`CaAlSi`$ and $`SrAlSi`$. The heat capacity data are useful to (i) measure the thermodynamic properties and confirm the bulk nature of superconductivity, (ii) compare the specific heat in the normal and superconducting states with the standard BCS theory,13 ; 8  and (iii) extract microscopic parameters, e.g. the superconducting gap. The application of pressure is of interest because an extraordinary large superconducting $`T_c`$ of 14 K was reported for a high pressure phase of the binary silicide $`CaSi_2`$.9

The heat capacity data clearly show the existence of bulk superconductivity in $`CaAlSi`$ and $`SrAlSi`$, consistent with previous magnetic measurements.1 ; 2  The temperature dependence of the specific heat of $`SrAlSi`$ fits well the BCS theory in the weak-coupling limit. Deviations from the BCS model in the case of $`CaAlSi`$ indicate a slightly stronger coupling. With the application of hydrostatic pressure $`T_c`$ increases in $`CaAlSi`$ but it decreases in $`SrAlSi`$. The opposite signs of the pressure coefficients of $`T_c`$ are explained in terms of the distinct difference between the electronic structures of both compounds.

## II Experimental

The pseudo-ternary intermetallic compounds, $`CaAlSi`$ and $`SrAlSi`$, have been synthesized by argon arc melting with appropriate amounts of $`Ca`$ (99 %), $`Sr`$ (99 %), $`Al`$ (99.99 %), and $`Si`$ (99.99 %) as described earlier.1  The x-ray powder diffraction spectra show the $`C32`$ hexagonal structure with lattice constants (a,c) = (4.189, 4.400 Å) for $`CaAlSi`$ and (a,c) = (4.220, 4.754 Å) for $`SrAlSi`$. Magnetization and transport measurements indicate the onset of superconductivity at 7.9 K and 5.0 K in $`CaAlSi`$ and $`SrAlSi`$, respectively.1 ; 4

The heat capacity was measured for both samples (typical mass of 25 mg) between room temperature and 2.2 K employing the Physical Property Measurement System (Quantum Design). The relative error of the measurement (i.e. the statistical error of the fitting procedure) was of the order of 0.1 to 0.2 %. The absolute accuracy of the measurement was tested in particular at low temperatures by comparing the measured heat capacity of high purity (99.9999 %) lead with the available data.10  The agreement of the lead measurement with the literature data was within the statistical data fluctuation.

The pressure effect on the superconducting transitions of $`CaAlSi`$ and $`SrAlSi`$ was investigated by ac susceptibility measurements at pressures up to 2 GPa. A dual coil system was mounted directly to the sample. The ac susceptibility was measured using the mutual inductance bridge (LR700, Linear Research). Pressure was generated in a beryllium-copper piston-cylinder clamp.11  The sample was mounted in a Teflon container filled with a 1:1 mixture of Fluorinert FC70 and FC77 as a hydrostatic pressure transmitting medium. The pressure was measured in situ at 7 K by monitoring the shift of the superconducting $`T_c`$ of a high purity (99.9999 %) lead manometer. The temperature above 45 K was measured by a thermocouple inside the Teflon container and, at low temperatures, by a germanium resistor built into the pressure cell close to the sample position. Data have been taken upon loading and unloading cycles.

## III Results and Discussion

### III.1 Heat capacity measurements

Figs. 1a and 1b show the low temperature heat capacity data, $`C_p/T`$ vs. $`T,`$ for $`CaAlSi`$ and $`SrAlSi`$, respectively. At zero magnetic field (open circles) $`C_p`$ shows a sharp peak at the superconducting transition as expected for bulk superconductors with a narrow transition width of 0.4 K ($`CaAlSi`$) and 0.2 K ($`SrAlSi`$). To compare with the normal state heat capacity the superconductivity was suppressed by a 7 Tesla magnetic field. The corresponding data are shown by open squares in Fig. 1. The high field data can be used to estimate the superconducting component of $`C_p`$ by subtracting the lattice contribution from the zero field $`C_p`$ below $`T_c`$. In general $`C_p`$ is given by the sum of the electronic ($`C_e`$) and the lattice ($`C_l`$) contributions,

$$C_p(T)=C_e(T)+C_l(T)$$
(1)
with $`C_e(T)=C_e^{(n)}(T)=\gamma T`$ in the normal state and $`C_e(T)=C_e^{(s)}(T)`$ in the superconducting state. The values of $`C_e^{(s)}(T)`$ in the BCS theory for weak coupling were numerically calculated and tabulated by Mühlschlegel.13  The characteristic coefficient $`\gamma `$ of $`C_e^{(n)}`$ is proportional to the density of states (DOS) at the Fermi energy, $`N(E_F)`$. The lattice heat capacity at low temperatures can be expanded in a power series of $`T`$ (Debye model) starting with the third order term, $`T^3`$. Normal state electronic and lattice contributions are commonly separated in the $`C_p/T`$ vs. $`T^2`$ plot where the coefficient $`\gamma `$ is obtained as the linear extrapolation of $`C_p/T`$ to zero temperature and the slope at low temperatures is used to calculate the Debye temperature, $`\mathrm{\Theta }_D`$. Using the high field data of $`C_p`$ we estimate $`\gamma =5.04mJ/molK,\mathrm{\Theta }_D=226K`$ and $`\gamma =5.48mJ/molK,\mathrm{\Theta }_D=218K`$ for $`CaAlSi`$ and $`SrAlSi`$, respectively. It is interesting to note that $`\gamma `$ of $`SrAlSi`$ is indeed larger than the value for $`CaAlSi`$ indicating that $`N(E_F)`$ increases with the substitution of $`Ca`$ by $`Sr`$ in qualitative agreement with the band structure calculations.7  However, the increase is moderate and far smaller than the factor of two predicted by the theory.

The electronic heat capacity in the superconducting state, $`C_e^{(s)}(T)`$, is compared with the BCS data (dotted line) in the insets of Fig. 1. The thermodynamic $`T_c`$ was estimated by an entropy conserving construction as $`T_c=7.7K`$ ($`CaAlSi`$) and $`T_c=4.9K`$ ($`SrAlSi`$). These values are slightly lower than the values determined from resistivity measurements or from the onset of the diamagnetic drop of the magnetic susceptibility. However, this is to be expected since the thermodynamic $`T_c`$ represents an averaged bulk value whereas resistivity as well as the onset of diamagnetism are both sensitive to the first percolating path of superconducting volume in the sample. The normalized electronic heat capacity, $`C_e^{(s)}(T)/C_e^{(n)}(T)`$, is displayed as a function of $`T/T_c`$. There is an excellent agreement with the BCS data in the case of $`SrAlSi`$ (Fig. 1b), i.e. the BCS model in weak coupling limit describes the superconducting state very well in the temperature range accessible in our experiments ($`T>2.2K`$). There are deviations of $`C_e^{(s)}(T)`$ from the BCS data in case of $`CaAlSi`$. The heat capacity falls below the BCS values at low $`T`$ but is enhanced close to $`T_c`$ resulting in a larger relative jump of $`C_p`$, $`\delta C_p(T_c)/\gamma T_c2`$ (BCS value = 1.43). These deviations may be explained by a stronger electron phonon coupling in $`CaAlSi`$ as compared to $`SrAlSi`$. The thermodynamic consistency of the heat capacity data was checked by integrating the entropy difference, $`(C_e^{(s)}(T)C_e^{(n)}(T))/T`$, from $`T=0`$ to $`T=T_c`$. For $`CaAlSi`$ this integral is close to zero as expected. For $`SrAlSi`$ the integration cannot be carried out with sufficient accuracy because of the experimental restrictions (the lowest accessible temperature, 2.2 K, is just little less than 50 % of $`T_c`$). The superconducting gap at zero temperature, $`2\mathrm{\Delta }`$, may be estimated from the heat capacity at low temperatures by plotting $`lnC_e^{(s)}(T)`$ vs. $`1/T`$ (Arrhenius plot). The low temperature slopes estimated from this plot yield values of $`2\mathrm{\Delta }=1.4meV`$ and $`2\mathrm{\Delta }=2.7meV`$ for $`SrAlSi`$ and $`CaAlSi`$, respectively. However, because of the limited temperature range of our data ($`T>`$2.2 K) these values are certainly underestimated and should be considered as a lower limit of the gap parameter. A more accurate estimate requires heat capacity measurements at far lower temperatures. The BCS value of the gap is related to $`T_c`$ by $`2\mathrm{\Delta }_{BCS}=3.53k_BT_c`$. For $`SrAlSi`$ with $`T_c=4.9K`$ we get $`2\mathrm{\Delta }_{BCS}=1.49meV`$ which is only slightly larger than our estimate. Taking into account the almost perfect agreement of $`C_e^{(s)}(T)`$ with the BCS function over the full experimental temperature range, the superconducting gap value for $`SrAlSi`$ can be assumed to be close to the BCS value of $`1.49meV`$. However, the calculated BCS gap parameter for $`CaAlSi`$, $`2\mathrm{\Delta }_{BCS}=2.34meV`$, is clearly smaller than the value of $`2.7meV`$ estimated as a lower limit of $`2\mathrm{\Delta }`$ from the Arrhenius plot of $`C_e^{(s)}(T)`$. This enhancement of the gap with respect to the BCS parameter is a further indication of an increased electron-phonon coupling. A summary of all relevant parameters and the comparison to the BCS weak coupling values is given in Table 1.

The temperature dependence of the heat capacity in the superconducting state provides strong evidence that the superconductivity found in $`SrAlSi`$ and $`CaAlSi`$ is well explained by the standard BCS theory for weak and enhanced coupling, respectively. Furthermore, the data confirm the bulk nature of the superconductivity in these pseudo-ternary $`C32`$ compounds.

It appears interesting to compare the present results with the superconductivity found at relatively high temperatures in binary compounds with similar structure. Superconductivity at 14 K was recently detected in $`CaSi_2`$ at pressures above 14 GPa.9  The (tetragonal) structure of this high pressure phase is very similar to the $`AlB_2`$ structure but the honeycomb planes of the $`Si`$ ions show still a nonvanishing buckling even at the highest pressures.14  No thermodynamic data such as heat capacity are known for this superconducting high pressure phase. In contrast, extensive research has been conducted to investigate the thermodynamic properties in the superconducting state of $`MgB_2`$. The heat capacity of $`MgB_2`$ shows an abnormal temperature dependence below $`T_c`$.15  Whereas the value of $`C_e^{(s)}(T)`$ at and right below $`T_c`$ is smaller than the corresponding BCS function it raises above the BCS data below about 20 K followed by an exponential drop at lower temperature. From this unusual behavior the existence of two superconducting gaps in $`MgB_2`$ was suggested and later confirmed by many alternative experiments. The present data for $`CaAlSi`$ and $`SrAlSi`$ show no evidence for a similar scenario within the experimental temperature range of $`T>`$2.2 K, i.e. $`T>0.44T_c`$ for $`SrAlSi`$ and $`T>0.27T_c`$ for $`CaAlSi`$. The existence of an anomaly at far lower temperature, similar to $`MgB_2`$, cannot be completely excluded from the current investigation. However, the almost perfect agreement of $`C_e^{(s)}(T)`$ with the BCS data in the case of $`SrAlSi`$ as well as the enhancement of $`C_e^{(s)}(T)`$ close to $`T_c`$ of $`CaAlSi`$ (instead of the depletion observed in $`MgB_2`$) suggest a less complex structure of the superconducting gap as compared to magnesium diboride.

### III.2 Pressure effect on the superconducting transitions in CaAlSi and SrAlSi

The effect of hydrostatic pressure on the superconducting $`T_c`$ is of special interest since it can help to reveal the intrinsic mechanisms of superconductivity. In $`MgB_2`$, for example, it could be shown that the negative pressure coefficient of $`T_c`$ was most compatible with the strong-coupling model (McMillen) of phonon mediated superconductivity so that alternative explanations appeared to be less favorable.12

The pressure effect on $`T_c`$ of $`CaAlSi`$ and $`SrAlSi`$ was investigated by measuring the real part of the ac susceptibility, $`\chi _{ac}^{}`$. The data in Fig. 2 show the shift of the diamagnetic signal of $`\chi _{ac}^{}`$ with increasing pressure. Surprisingly, the pressure effect on $`CaAlSi`$ and $`SrAlSi`$ is opposite in sign. The superconducting transition temperature (measured at the onset of the diamagnetic signal) of $`CaAlSi`$ increases with pressure whereas $`T_c`$ of $`SrAlSi`$ clearly decreases. The $`T_c`$’s as a function of pressure for both compounds are shown in Fig. 3. The increase of $`T_c`$ in $`CaAlSi`$ (Fig. 3a) is slightly nonlinear with an initial slope of $`0.21K/GPa`$. In contrast, Fig. 3b shows an almost perfect linear decrease of $`T_c`$ for $`SrAlSi`$ with a coefficient of $`0.12K/GPa`$. Although the magnitude of the relative pressure coefficient, $`dln(T_c)/dp`$, is almost the same for both compounds, the opposite sign deserves further consideration.

For a qualitative discussion we use the BCS equation for $`T_c`$8

$$k_BT_c=1.13\text{h}\text{ }\text{ }\omega _D\mathrm{exp}\left(1/\lambda \right)$$
(2)
h  $`\omega _D`$ is a measure of the characteristic energy scale for the phonons and $`\lambda =V_0N(E_F)`$. $`V_0`$ is the matrix element of the effective interaction. Assuming that the pressure dependence of $`V_0`$ is negligible, the pressure coefficient of $`T_c`$ is derived as

$$\frac{d\mathrm{ln}T_c}{dp}=\frac{d\mathrm{ln}\omega _D}{dp}+\frac{1}{V_0N(E_F)}\frac{d\mathrm{ln}\left[V_0N(E_F)\right]}{dp}$$
(3)
There are obviously two contributions to the pressure coefficient of $`T_c`$: one from the sole phonon system, $`\frac{d\mathrm{ln}\omega _D}{dp}`$, and the second one involving the electronic system. In the weak-coupling limit of the BCS model both terms contribute to the total pressure coefficient. The effect of pressure on the phonon system usually results in an increase of the average phonon energy, i.e. in a “hardening” of phonon modes. Only in rare cases, close to a structural transition, phonon “softening” may occur but usually only for some particular modes. Therefore, the first term in equation (3) yields a positive contribution to the total pressure coefficient. The pressure effect on the density of states, however, is more complex. The compression of the lattice by pressure in general causes an increase of the bandwidth which in turn will result in an average decrease of $`N(E)`$ since the total number of states in a band is fixed. In addition, $`V_0`$ tends to decrease with pressure since it is inversely proportional to the average square of the phonon frequency. This negative contribution to $`\frac{d\mathrm{ln}T_c}{dp}`$ will compete with the phonon hardening effect. However, the application of pressure may also cause slight changes in the band structure and/or a shift of the Fermi energy, $`E_F`$, and a change of $`N(E_F)`$. Therefore, within the BCS theory the pressure coefficient of $`T_c`$ can be of either sign, depending on which of the contributions to equation (3), positive or negative, dominate.

In both compounds, $`CaAlSi`$ and $`SrAlSi`$, the phonon contribution to the pressure coefficient is expected to be positive. The magnitude of $`\frac{d\mathrm{ln}\omega _D}{dp}`$ is not known so far but it can be estimated from Raman or infrared spectroscopy at high pressures. The electronic term in equation (3) needs a careful consideration. Recent band structure calculations have shown striking differences in $`N(E)`$ in $`CaAlSi`$ and $`SrAlSi`$.7  In $`CaAlSi`$ the Fermi energy lies in a region where $`N(E)`$ is flat and relatively insensitive to small variations of $`E_F`$. This results in a small electronic contribution to the total pressure coefficient of $`T_c`$. In $`CaAlSi`$ the positive phonon term obviously outweighs the electronic term which results in the observed positive $`\frac{dT_c}{dp}`$ as shown in Fig. 3a. The situation is different, however, in $`SrAlSi`$. According to the calculations7  the density of states forms a very narrow peak close to the Fermi energy. This makes $`N(E_F)`$ extremely susceptible to external perturbations. Any small change induced by pressure (band broadening, shift of $`E_F`$) may result in a stronger decrease of $`N(E_F)`$. The negative electronic contribution to the pressure coefficient of $`T_c`$ therefore can dominate over the lattice term and the total $`\frac{dT_c}{dp}`$ becomes negative as observed in our experiments (Fig. 3b).

Although the above discussion is qualitative in nature it indicates that the opposite pressure coefficients of $`T_c`$ in $`CaAlSi`$ and $`SrAlSi`$ uniquely reflect the peculiarities of their electronic structure. They provide an indirect proof of the distinct differences in the density of states at $`E_F`$ in both compounds, $`CaAlSi`$ and $`SrAlSi`$, as suggested by band structure calculations. Furthermore, the positive pressure coefficient observed in $`CaAlSi`$ cannot be explained by a sole change of the density of states induced by pressure but the contribution due to phonon hardening has to be taken into account. For a more quantitative comparison additional experimental as well as theoretical work is needed. The measurement of phonon frequencies as a function of pressure should provide a better estimate of the phonon effect on $`\frac{d\mathrm{ln}T_c}{dp}`$. Improved band structure calculations involving the effects of lattice compression and including the phonon system (see for example ref. 16) may help to get a more adequate understanding of the superconductivity in the pseudo-binary C32 intermetallic compounds.

## IV Summary and Conclusions

We have investigated the superconductivity in the $`C32`$ intermetallic compounds $`CaAlSi`$ and $`SrAlSi`$ by heat capacity and high-pressure ac susceptibility measurements. Several parameters characterizing the normal and superconducting states are extracted from the heat capacity. The coefficients of the normal state electronic specific heat are estimated as $`\gamma =5.04mJ/molK^2`$ and $`\gamma =5.42mJ/molK^2`$ for $`CaAlSi`$ and $`SrAlSi`$, respectively. The values of the superconducting gap are calculated from the electronic heat capacity in the superconducting state as $`2.7meV`$ ($`CaAlSi`$) and $`1.4meV`$ ($`SrAlSi`$). The parameters are consistent with the predictions of the BCS theory in the weak-coupling limit for $`SrAlSi`$ and for enhanced coupling in the case of $`CaAlSi`$. The bulk nature of superconductivity in $`CaAlSi`$ and $`SrAlSi`$ is confirmed.

The effect of hydrostatic pressure on $`T_c`$ is positive in the case of $`CaAlSi`$ ($`\frac{d\mathrm{ln}T_c}{dp}=0.026GPa^1`$) but negative for $`SrAlSi`$ ($`\frac{d\mathrm{ln}T_c}{dp}=0.024GPa^1`$). This opposite tendency of the pressure shift of $`T_c`$ reflects the distinct differences in the density of states near the Fermi energy in both compounds. The high pressure data also underline the importance of the phonon frequency (or ionic mass) for understanding the superconductivity in the system $`AAlSi`$ and, in particular, the positive pressure effect observed in $`CaAlSi`$. With the increasing mass of the alkaline-earth ion $`T_c`$ is systematically suppressed, from 7.7 K ($`Ca`$) to 4.9 K ($`Sr`$) and to non-superconducting above 2 K ($`Ba`$). This tendency cannot be explained as a sole density-of–states effect since our data show that $`N(E_F)`$ increases when $`Ca`$ is replaced by $`Sr`$, in qualitative agreement with the band structure calculations.

###### Acknowledgements.

This work is supported in part by NSF Grant No. DMR-9804325, the T.L.L. Temple Foundation, the John J. and Rebecca Moores Endowment, and the State of Texas through the TCSAM at the University of Houston and at Lawrence Berkeley Laboratory by the Director, Office of Energy Research, Office of Basic Energy Sciences, Division of Materials Sciences of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.