text:
# Structure Space of Model Proteins: A Principal Component Analysis

## I INTRODUCTION

Proteins fold into specific structures to perform their biological function. Despite the huge diversity in their functions, evolutional paths, structural details, and sequences, the vast majority of proteins adopt only a small number ($`1000`$) of folds (“topology”).  This observation has intrigued a number of authors and lead to the concept of designability.  The designability of a structure is defined to be the number of sequences that have that structure as their unique lowest-energy state.  It has been shown in various model studies that structures differ drastically in their designability; a small number of highly designable structures emerge with their associated number of sequences much larger than the average.  Highly designable structures are also found to possess other protein-like properties, such as thermodynamic stability,  fast folding kinetics,  and tertiary symmetry.  These results suggest that there may be a designability principal behind nature’s selection of protein folds; these small number of folds were selected because they are readily designed, stable against mutations, and thermodynamically stable.

Why are some structures more designable than others? How do we identify highly designable structures? Finkelstein and co-workers argued that certain motifs are easier to stabilize and thus more common because they either have lower (e.g. bending) energies or have unusual energy spectra over random sequences.  Govindarajan and Goldstein suggested that the topology of a protein structure should be such that it is kinetically “foldable”.  More recently, it was noted that an important clue resides in the distribution of structures in a suitably defined structure space, with highly designable structures located in regions of low density.  In particular, within a hydrophobic model, Li et al. showed that the distribution of structures is very nonuniform, and that the highly designable structures are those that are far away from other structures.  However, identifying highly designable structures still remains a tedious task, requiring either full enumeration or sufficiently large sampling of both the structure and the sequence spaces, making studies of large systems prohibitive.

In this paper, we investigate the properties of the structure space of the hydrophobic model of Li et al., starting from a Principal Component Analysis (PCA). We show that while the distribution of the structures is not uniform, it can be approximated as a cloud of points centered on a single peak. The principal directions of this cloud are almost coincident with those obtained by rotation into Fourier space; the coincidence is in fact exact for the subset of cyclic structures. An interesting feature is that the eigenvalues of PCA, describing the extent of the density cloud along the principal axis, vary continuously with the Fourier label $`q`$, with a minimum at $`q=\pi `$ corresponding to alternating patterns. The continuity of the eigenvalues suggests an expansion around $`q=\pi `$, which leads to an analytical conjecture for the density of structures in the $`N`$-dimensional binary space. Assuming the validity of this conjecture in more general models, it provides a means of estimating density, and hence indirectly designability, of structures by simply analyzing their sequences, without the need for extensive enumerations of other possible structures.

The rest of the paper is organized as follows. In Section II we review the hydrophobic model and the designabilities of structures. In Section III we discuss the methods and the results of PCA applied to the structure space, and relate the density and designability of a structure to its projections onto the principal axes. In Section IV we demonstrate that Fourier transformation provides a very good approximation to PCA, and show that in fact the two procedures are equivalent for the subset of cyclic structures. As a comparison with real structures, in Section V we introduce and study an ensemble of pseudo-structures constructed by a Markovian process. Finally, in Section VI we synthesize the numerical results of PCA analysis, and develop a conjecture for the density of points in structure space.

## II THE HYDROPHOBIC MODEL

We start with a brief review of the hydrophobic model of Li et al.  and the designabilities of structures. Model sequences are composed of two types of amino acids, H and P. Each sequence $`\{h_i\}`$ (for $`i=1,2,\mathrm{},N`$) is represented by a binary string or vector, with $`h_i=0`$ for a P-mer and $`h_i=1`$ for an H-mer. We take the polymer length $`N=36`$, for which there are $`2^{36}`$ sequences. Each of these sequences can fold into any one of the many compact structures on a $`6\times 6`$ square lattice (Fig. 1). There are $`57,337`$ such compact structures unrelated by rotation and mirror symmetries. In the hydrophobic model, the only contribution to the energy for a sequence folded into a structure is the burial of the H-mers in the sequence into the core of the structure. So if one represents a structure by a binary string or vector, $`\{s_i\}`$, for $`i=1,2,\mathrm{},36`$, with $`s_i=0`$ for the surface sites and $`s_i=1`$ for the core sites (Fig. 1), the energy is

$$E=\underset{i=1}{\overset{N}{}}h_is_i,$$
(1)
where $`h_i`$ is the sequence vector.

The designability of a structure is defined as the number of sequences that have the structure as their unique lowest-energy state. To obtain an estimate for designabilities of structures, we randomly sampled $`50,000,000`$ sequences and found the unique lowest-energy structure, if any, for each of them. In Fig. 2, we plot the histogram of designabilities, i.e. number of structures with a given designability. Note that we have normalized designability so that its maximum value of $`2981`$ is scaled to one. In this paper, we define highly designable structures to be the top one percent of designable structures (structures with nonzero designability), which means $`307`$ structures with a designability larger than $`0.47`$.

In the hydrophobic model, both sequences and structures can be regarded as points in a 36-dimensional binary space, or corners of a hypercube in a Euclidean space of similiar dimension. In this representation, the lowest-energy state of a sequence is simply its nearest structure point.  Designabilities can then be obtained by constructing Voronoi polyhedra around all points corresponding to structures in this space; the designability of each structure is then the number of sequence points that fall within the corresponding Voronoi polytope (Fig. 3). Structures in the lower density regions have larger Voronoi polytopes and higher designability. Understanding how the structure points are distributed in this 36-dimensional space can thus help us address questions concerning designability. In the next section we examine the distribution of the structure points via the method of PCA.

## III PRINCIPAL COMPONENT ANALYSIS

First, let us note that while sequences are uniformly distributed in the 36-dimensional hypercube, structures are distributed on a 34-dimensional hyperplane because of the following two geometrical constraints. The first constraint on structure vectors comes from the fact that all compact structures have the same number of core sites, and thus

$$\underset{i=1}{\overset{36}{}}s_i=16.$$
(2)
The second constraint is that since the square lattice is bipartite, and any compact structure traverses an equal number of ‘black’ and ‘white’ points,

$$\underset{i=1}{\overset{36}{}}(1)^is_i=0.$$
(3)
Next, let us define the covariance matrix of the structure space as

$$C_{i,j}=s_is_js_is_j,$$
(4)
where $`i,j=1,2,\mathrm{},36`$, and the average is over all the $`57,337`$ possible $`(s_1,s_2,\mathrm{},s_{36})`$ for compact structures. The $`36\times 36`$ covariance matrix is symmetric $`C_{i,j}=C_{j,i}`$, and also satisfies the condition $`C_{i,j}=C_{37i,37j}`$. The latter is due to the reverse-labeling degeneracy of the structure ensemble, since if the string $`(s_1,s_2,\mathrm{},s_{36})`$ is in this ensemble, then its reverse $`(s_{36},s_{35},\mathrm{},s_1)`$ is also included. This symmetry implies that if $`(v_1,v_2,\mathrm{},v_{36})`$ is an eigenvector of the matrix $`C_{i,j}`$, then $`(v_{36},v_{35},\mathrm{},v_1)`$ is also an eigenvector with the same eigenvalue. Therefore, for every eigenvector of $`C_{i,j}`$ we have either $`v_j=v_{37j}`$ or $`v_j=v_{37j}`$.

As depicted in Fig. 4, the matrix $`C_{ij}`$ is peaked along the diagonal and decays off-diagonally with short range correlations. This feature reflects a general property of compact self-avoiding walks; if a monomer is in the core (on the surface) the neighboring monomers along the chain have enhanced probability to be in the core (on the surface). Another characteristic of $`C_{ij}`$ is that it is almost a function of $`|ij|`$ only, i.e. $`C_{ij}F(|ij|)`$, barring some small end and parity effects. We expect this feature of approximate translational invariance to be generic beyond the $`6\times 6`$ lattice model studied here. We also looked at the covariance matrix for the subset of highly designable structures. While qualitatively similar, it tends to decay faster off-diagonally than that of all structures. This is attributed to the fact that highly designable structures tend to have more frequent transitions between core and surface sites.

For PCA of structure space, the matrix $`C_{ij}`$ is diagonalized to obtain its eigenvectors $`\{\stackrel{}{v}^{(k)}\}`$, and the corresponding eigenvalues $`\{\lambda _k\}`$ for $`k=1,2,\mathrm{},36`$, which are shown in Fig. 5. The two zero eigenvalues ($`\lambda _1=\lambda _2=0`$) result from the constraints in Eqs. (2) and (3), with the corresponding eigenvectors of $`v_i^{(1)}=1`$, and $`v_i^{(2)}=(1)^i`$ for $`i=1,2,\mathrm{},36`$, respectively. The remaining 34 nonzero eigenvalues range smoothly from zero to one, making any further dimensional reduction not obvious. For comparison, the 36 eigenvalues of the uniformly distributed points of sequence space are all the same ($`\lambda =1/4`$). (It is easy to show that the covariance matrix for the sequence space is $`C_{ij}=\delta _{ij}/4`$.) On the other hand, a uniform distribution on the 34-dimensional hyperplane where the structure points reside would result in 34 identical eigenvalues of $`360/13770.26`$.

Identification of the principal axes and eigenvalues does not necessarily provide information about the distribution of points in space. To examine the latter, we first project each structure vector onto its components along the eigenvectors. Using the rotation matrix $`R_{ki}`$ that diagonalizes the covariance matrix, the component $`y_k`$ of the structure vector along principal axis $`k`$ is obtained as

$$y_k=\underset{i=1}{\overset{36}{}}\left(s_is_i\right)R_{ki}.$$
(5)
Interestingly, we find that along each of the principal directions, the distribution of components is a bell-shaped function with a single peak at zero. Such distributions can then be well approximated by Gaussians whose variances are the corresponding eigenvalues $`\lambda _k`$, i.e.

$$\rho _k(y_k)\frac{1}{\sqrt{2\pi \lambda _k}}e^{\frac{y_k^2}{2\lambda _k}}.$$
(6)
In Fig. 6 we show the distribution of projections $`y_k`$ on two principal axes $`k=16`$ and $`k=36`$, along with the corresponding Gaussian distributions.

Equation (6) provides a good characterization of the density of structures in the $`N`$ dimensional space. Highly designable structures are expected to lie in regions of this space where the density of structures is small, while the number of available sequences is large. Let us consider a structure characterized by a vector $`\stackrel{}{y}`$. If the density of structural points in the vicinity of this point is $`\rho _{str}(\stackrel{}{y})`$, the number of available structures in a volume $`V`$ around this point is $`V\rho _{str}(\stackrel{}{y})`$. Neglecting various artifacts of discreteness, the volume of the Voronoi polyhedron (see Fig. 3) around this point is given by $`V(\stackrel{}{y})1/\rho _{str}(\stackrel{}{y})`$. The designability is the number of structures within this volume, and estimated as $`\rho _{seq}(\stackrel{}{y})/\rho _{str}(\stackrel{}{y})`$, where $`\rho _{seq}(\stackrel{}{y})`$ is the density of sequences. The sequence density is in fact uniform in the $`N`$-dimensional space. The structure density can be approximated as the product of Guassians along the principal projections, and thus

$$\mathrm{Designability}\frac{\rho _{seq}(\stackrel{}{y})}{\rho _{str}\left(\stackrel{}{y}\right)}\underset{k=3}{\overset{36}{}}\frac{1}{\rho _k\left(y_k\right)}\mathrm{exp}\left[\underset{k=3}{\overset{36}{}}\frac{y_k^2}{2\lambda _k}\right](\stackrel{}{y}).$$
(7)
We have neglected various proportionality constants in the above equation, leading to the quantity $`(\stackrel{}{y})`$ which is our estimator for designability. In Fig. 7, the estimate $``$ is plotted against the actual designability for all designable structures. There is a reasonably good, but by no means perfect, correlation between the designability and the estimator $``$. The structures with the top one percent value of $``$ include 39% of the highly designable structures.

## IV FOURIER DECOMPOSITION AND CYCLIC STRUCTURES

In discussing Fig. 4, we already noted that the covariance matrix $`C_{ij}`$ is approximately a function of $`|ij|`$, with corrections due to end effects. If this were an exact symmetry, the matrix would be diagonal in the Fourier basis. Even in the presence of the end effects, Fourier decomposition provides a very good approximation to the eigenvectors and eigenvalues of PCA, as demonstrated below. For each structure vector $`\{s_j\}`$, the Fourier components are obtained from

$$S_q=\frac{1}{\sqrt{N}}\underset{j=1}{\overset{N}{}}e^{iqj}\left(s_js_j\right),$$
(8)
where $`q=2\pi \alpha /N`$, with $`\alpha =0,1,\mathrm{},N1`$. The average value of $`s_j`$ is subtracted for convenience. With this subtraction, the two constraints in Eqs. (2) and (3) correspond to two zero modes in Fourier space, as $`S_0=0`$ and $`S_\pi =0`$, and since $`\{s_j\}`$ are real $`S_q^{}=S_q`$.

The covariance matrix in the Fourier space is

$$S_qS_q^{}^{}=\frac{1}{N}\underset{j,j^{}=1}{\overset{N}{}}e^{i(qjq^{}j^{})}C_{jj^{}},$$
(9)
and is both real and symmetric (since $`C_{jj^{}}=C_{j^{}j}`$). If $`C_{jj^{}}`$ is translationally invariant, i.e. $`C_{jj^{}}=F(|jj^{}|modN)`$, Eq. (9) becomes

$$S_qS_q^{}^{}=\delta _{q,q^{}}\lambda _q,$$
(10)
where

$$\lambda _q=\underset{k=0}{\overset{N1}{}}e^{iqk}F(k)=\left|S_q\right|^2,$$
(11)
are the diagonal elements of the diagonalized matrix in Eq. (10), and hence the eigenvalues of $`C_{jj^{}}`$. Note that because the matrix is real–symmetric, its eigenvalues appear in pairs, i.e.

$$\lambda _q=\lambda _q.$$
(12)
Since our covariance matrix is not fully translationally invariant, $`S_qS_q^{}^{}`$ is not diagonal. However, as shown in Fig. 8a, its off diagonal elements are very small. As required by symmetry, the diagonal elements form pairs of identical values. These diagonal elements, plotted versus the index $`\alpha `$ in Fig. 8b, should provide a good approximation to the eigenvalues obtained by PCA. This is corroborated in Fig. 9a, where we compare $`|S_q|^2`$ with the true eigenvalues of the covariance matrix $`C_{jj^{}}`$.

Finally, we note that the end effects that mar the translational invariance of the covariance matrix are absent in the subset of cyclic structures. Any structure whose two ends are neighboring points on the lattice can be made cyclic by adding the missing bond. Any one of the $`N=36`$ bonds on the resulting closed loop can be broken to generate an element of the original set of structures, and the corresponding structure strings are cyclic permutations of each other. Thus, the covariance matrix $`C_{cyclic}(j,j^{})`$ of the set of all cyclic structures is translationally invariant. In our model of $`6\times 6`$ compact polymers, there are a total of $`36\times 276`$ cyclic structures. The Fourier transform of their covaraince matrix is diagonal as expected, with diagonal elements depicted in Fig. 9b. The corresponding Fourier eigenvalues are quite close to the eigenvalues of the full matrix obtained in the PCA (Fig. 9b). Thus the end effects do not significantly modify the correlations, and this is specially true for the smallest eigenvalues which make the largest contributions to the density in Eq. 7.

## V A MARKOVIAN ENSEMBLE of PSEUDO–STRUCTURES

The geometry of the lattice, and the requirement of compactness constrain the allowed structure strings of zeros and ones in a non-trivial fashion. In our estimation of designabilities so far, we have focused on the covariance matrix which carries information only about two point correlations along these strings. In principle, higher point correlations may also be important, and we may ask to what extent the covariance matrix contains the information about the structures’ designabilities? As a preliminary test, we performed a comparative study with an artificial set of strings, not corresponding to real structures, but constructed to have a covariance matrix similar to true structures on the $`6\times 6`$ lattice.

Specifically, we generated a set of random strings $`\{\stackrel{}{t}\}`$, of zeros and ones of length 36, using a third order Markov process as follows. For each string, the first element $`t_1`$ is generated with probability $`P(t_1=1)=s_1`$, where $`s_1`$ is the fraction of the true structure strings with $`s_1=1`$. The second element $`t_2`$ is generated according to a transition probability $`P(t_1t_2)`$ which is taken to be the “conditional probability” $`P(s_2|s_1)`$ extracted from the true structure strings. The third point $`t_3`$ is generated according to a transition probability $`P(t_1t_2t_3)`$ which is the “conditional probability” $`P(s_3|s_1s_2)`$ extracted from the true structure strings. All the remaining points $`t_j`$, $`j=4,5,\mathrm{},36`$, are generated according to the transition probabilities $`P(t_{j3}t_{j2}t_{j1}t_j)`$ equal to the true “conditional probabilities” $`P(s_j|s_{j3}s_{j2}s_{j1})`$ of actual structures. Sequences that do not satisfy the global constrains of Eqs. (2) and (3) are thrown out. For every Markov string generated, we also put its reverse in the pool, unless the string is its own reverse.

The above Markovian ensemble has a covariance matrix, and corresponding eigenvalues, very similar to those of the true structures, as shown in Fig. 10. We then calculated the designabilities of these “pseudo-structures” using Eq. (1) by uniformly sampling $`5\times 10^7`$ random binary sequences. The histogram of the designabilities (Fig. 11) is qualitatively similar to that of the true structures (Fig. 2). Next we constructed the designability estimator $``$ (Eq. (7)) for the pseudo-structures, using the eigenvalues and eigenvectors of their covariance matrix. The quantity $``$ is plotted versus designability in Fig. 12 for all the artificial pseudo-structures with non-zero designability. The pseudo-structures with the top one percent value of $``$ include 60% of the highly designable psuedo-structures.

These results suggest that a considerable amount of information about the designability is indeed contained in the two point correlations of the string. The designability estimator, Eq. (7) in fact does a somewhat better job in the case of pseudo-structures generated according to short-range Markov rules.

## VI CONCLUSIONS

One of the most intriguing properties of compact structures, which emerged from early extensive enumeration studies, is that designabilites range over quite a broad distribution of values. Such a large variation in designability is a consequence of a non-uniform distribution of structure vectors, with highly designable structures typically found in regions of low density. However, our study of $`6\times 6`$ lattice structures using PCA indicates that the non-uniform density actually has a rather simple form that can be well approximated by a single multi-variable Gaussian, as in Eq. (7). Since this method of estimating structure designability is based only on the overall distribution of structures, it can be a useful tool in cases where there is not enough computational power to enumerate the whole structure space and calculate the designability. To obtain an accurate enough covariance matrix requires only a uniform sampling of the structure space.

We can also attempt to use the numerical results as a stepping stone to a more analytical approach for calculating the density of structures. First, we note that the covariance matrix for all structures is rather similar to that of the subset of cyclic structures, and that for the latter PCA is equivalent to Fourier decomposition. Second, we observe that the multi-variable Gaussian approximation to structure density in Eq. (7) is most sensitive to the eigenvalues that are close to zero. In terms of Fourier components, these are eigenvalues corresponding to values of $`q`$ close to $`\pi `$, and related to the constraint of Eq. (3). There is also a zero eigenvalue for $`q=0`$, related to condition (2). However, the latter global constraint appears not to have any local counterpart, as there is a discontinuity in the eigenvalues close to zero. Third, the continuity of the eigenvalues as $`q\pi `$, along with the symmetry of Eq. (12), suggests as expansion of the form $`\lambda _q=K(q\pi )^2+𝒪\left((q\pi )^4\right)`$. Indeed the numerical results indicate that the important (smaller) eigenvalues can well be approximated by $`K(q\pi )^2`$, with $`K0.04`$.

With this approximation, the designability estimate of Eq. (7) becomes

$$\left(\left\{\stackrel{}{s}\right\}\right)\mathrm{exp}\left[\underset{q}{}\frac{\left|S_q\right|^2}{2K(q\pi )^2}\right]=\mathrm{exp}\left[\frac{1}{2K}\underset{i,j=1}{\overset{N}{}}(1)^is_iJ_N\left(|ij|\right)(1)^js_j\right].$$
(13)
The first form in the above equation expresses the estimate in terms of the Fourier modes of the structure string, while the second term is directly in terms of the elements $`\{s_i\}`$. The function $`J_N\left(|ij|\right)`$ is the discrete Fourier transform of $`1/q^2`$, which for large $`N`$ behaves as $`|ij|`$. Equation (13) is thus equivalent to the Boltzmann weight of a set of unit charges on a discrete line of $`N`$ points marked by parity. The charges on the sublattice of the same parity attract each other with a potential $`J_N(r)`$, while those on different sublattices repel. Such an interaction gives a larger weight (and hence designability) to configurations in which the charges alternate between the core and surface sites, as observed empirically.

In would be revealing to see how much of the above results, developed on the basis of a lattice hydrophobic model, can be applied to real protein structures. One could use the exposure level of residues to the solvent in building up the structure vectors. Current methods deal with structure strings of a fixed length, equal to the dimension of the structure space. Since real proteins have different lengths, there is a need for a scaling method to handle them all together. Our study shows that the two point correlations of structure vectors are approximately translationally invariant, and can be captured by Fourier analysis. This suggests the possibility of casting the density of points in structure space in universal functional forms dependent only of a few parameters encoding the properties of the underlying polymers. If so, it would be possible to provide good estimates for designability with polymers of varying length, without the need for extensive numerical computations.

## VII ACKNOWLEDGEMENT

We thank G. Gatz, N. Butchler and E. Domany for useful discussions. The work was initiated, and partly completed, during the Program on Statistical Physics and Biological Information at the Institute for Theoretical Physics, UCSB, supported in part by the NSF Contract No. PHY99-07949. We acknowledge support by the NSF Grant No. DMR-01018213 (MY and MK), and an ITP Graduate Fellowship (MY).

FIGURE CAPTIONS

1. A possible compact structure on the $`6\times 6`$ square lattice. The 16 sites in the core region, enclosed by the dashed lines, are indicated by 1’s; the 20 sites on the surface are labeled by 0’s. Hence this structure is represented by the string 001100110000110000110011000011111100. Note that each ‘undirected’ open geometrical structure can be represented by two ‘directed’ strings, starting from its two possible ends (except for structures with reverse-labeling symmetry where the two strings are identical). It is also possible for the same string to represent different structures which are folded differently in the core region. For the $`6\times 6`$ lattice of this study, there are $`26929`$ such ‘degenerate’ structures, which are by definition non-designable.
2. Number of structures with a given designability versus relative designability for the $`6\times 6`$ hydrophobic model. The data is generated by uniformly sampling $`5\times 10^7`$ strings from the sequence space. The designability of each structure is normalized by the maximum possible designability.
3. Schematic representation of the 36-dimensional space in which sequences and structures are vectors or points. Sequences, represented by dots, are uniformly distributed in this space. Structures, represented by circles, occupy only a sparse subset of the binary points and are distributed non-uniformly. The sequences lying closer to a particular structure than to any other, have that structure as their unique ground state. The designability of a structure is therefore the number of sequences lying entirely within the Voronoi polytope about that structure.
4. Covariance matrix $`C_{ij}`$ of all compact structures of the $`6\times 6`$ square.
5. Eigenvalues of the covariance matrix for the structure vectors (circles), and for all points in sequence space (crosses).
6. Distributions of projections $`y_k`$ onto principal axes $`k=16`$ (a), and $`k=36`$ (b), for all 57337 structures (dots). Also plotted are Gaussian forms with variances $`\lambda _{16}`$ and $`\lambda _{36}`$, respectively (dashed lines).
7. The estimate $``$ (Eq. (7)) versus scaled designability for all designable structures on the $`6\times 6`$ square.
8. (a) The Fourier transformed covariance matrix $`S_qS_q^{}^{}`$ (Eq. (9)); and (b) its diagonal elements $`|S_q|^2`$.
9. (a)The diagonal elements $`|S_q|^2`$ (dots) plotted together with the eigenvalues of the covariance matrix (pluses). (b) Diagonal elements of the Fourier transformed $`C_{cyclic}`$, which are also the eigenvalues of the covariance matrix of cyclic structures (pluses). Eigenvalues of the covariance matrix for all structures are indicated by stars.
10. Eigenvalues of the covariance matrix for the structures generated by the Markov model (circles), and that of the true structure space (pluses).
11. Number of pseudo-structures with a given designability versus designability for the pseudo-structure strings randomly generated using the Markov model. The data is generated by uniformly sampling $`5\times 10^7`$ binary sequence strings. The designability of each pseudo-structure is normalized by the maximum possible designability.
12. The quantity $``$ versus designability for all designable pseudo-structures generated by the Markov model.


format: text

text:
# 1 Introduction

## 1 Introduction

The radiative decay $`\overline{B}X_s\gamma `$ is known to be extremely sensitive to the structure of fundamental interactions at the electroweak scale. It is dominantly generated by the Flavour Changing Neutral Current (FCNC) decay $`bs\gamma `$ that does not arise at the tree level in the Standard Model (SM). The leading order SM diagrams are shown in fig. 1.

Many possible non-standard contributions (e.g., SUSY one-loop diagrams) are of the same order in electroweak interactions. They might remain important even for relatively heavy exotic particles. Consequently, $`bs\gamma `$ imposes severe constraints on extensions of the SM (see, for instance, ).

The inclusive branching ratio BR$`[\overline{B}X_s\gamma ]`$ has been measured so far by CLEO , BELLE  and ALEPH . The most accurate result is the one of CLEO, where photons with energies down to 2.0 GeV are included. Extrapolation towards lower photon energy cutoffs is performed following the phenomenological models of refs. .

When the photon energy cutoff is chosen to be 1.6 GeV in the $`\overline{B}`$-meson rest frame, the experimental world average reads111 By convention, contributions to $`\overline{B}X_s\gamma `$ from intermediate real $`\psi `$ and $`\psi ^{}`$ are treated as background, while all the other $`c\overline{c}`$ states are included.

$$\mathrm{BR}[\overline{B}X_s\gamma (E_\gamma >1.6\mathrm{GeV})]_{\mathrm{exp}}=(3.12\pm 0.41)\times 10^4.$$
(1.1)
Within $`1\sigma `$, it matches the SM prediction

$$\mathrm{BR}[\overline{B}X_s\gamma (E_\gamma >1.6\mathrm{GeV})]_{\mathrm{SM}}=(3.57\pm 0.30)\times 10^4.$$
(1.2)
One can see that the experimental and theoretical uncertainties are close in size. Without the inclusion of the Next-to-Leading Order (NLO) QCD corrections, the theoretical uncertainty in eq. (1.2) would be around three times larger, and the constraints on new physics — much weaker.

The program of the NLO calculation was outlined by Stefan Pokorski, Manfred Münz and us in the article . At that time, the only known results were the Leading Order (LO) ones that suffered from large scale uncertainties . We analyzed these uncertainties in detail, and enumerated calculations that still had to be done in the NLO case. Very recently, the last element of this NLO program has been completed . In parallel to the QCD calculations, progress was being made in evaluation of the electroweak corrections, non-perturbative effects, as well as in collecting and analyzing the experimental data.

In the present paper, we summarize all the contributions to the NLO QCD calculation of BR$`[\overline{B}X_s\gamma ]`$, and discuss perspectives for further improvement of the theoretical accuracy. In particular, we point out the interplay between charm-quark mass uncertainties in the perturbative calculation and non-perturbative effects.

Our article is organized as follows. The next section is devoted to a brief description of the history of perturbative calculations of QCD effects in $`bs\gamma `$. In section 3, we summarize the electroweak corrections. Non-perturbative effects are discussed in section 4. The main theoretical uncertainties and possibilities for their elimination are the subject of section 5. Section 6 contains our conclusions.

## 2 The LO and NLO QCD calculations

In a certain range of photon energy cutoffs, the width of the hadronic decay $`\overline{B}X_s\gamma `$ is well approximated by the perturbative decay width

$$\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]=\mathrm{\Gamma }[bs\gamma ]+\mathrm{\Gamma }[bs\gamma g]+\mathrm{}.$$
(2.1)
Arguments that support such a statement will be discussed in section 4. Until then, we shall restrict our discussion to the perturbative quantity (2.1).

The framework for all the renormalization-group-improved perturbative analyses of $`bX_s^{\mathrm{parton}}\gamma `$ is set by the effective Lagrangian

$$_{\mathrm{eff}}=_{\mathrm{QCD}\times \mathrm{QED}}(u,d,s,c,b)+\frac{4G_F}{\sqrt{2}}V_{ts}^{}V_{tb}\underset{i=1}{\overset{8}{}}C_i(\mu )Q_i+\mathrm{}.$$
(2.2)
It is obtained from the underlying theory (SM in our case) by decoupling of all the particles that are much heavier than the $`b`$-quark. The Wilson coefficients $`C_i(\mu )`$ play the role of coupling constants at the vertices $`Q_i`$. The generic structure of the operators $`Q_i`$ is as follows:

$$Q_i=\{\begin{array}{cc}(\overline{s}\mathrm{\Gamma }_ic)(\overline{c}\mathrm{\Gamma }_{\underset{¯}{i}}^{}b),\hfill & i=1,2,\hfill \\ (\overline{s}\mathrm{\Gamma }_ib)_q(\overline{q}\mathrm{\Gamma }_{\underset{¯}{i}}^{}q),\hfill & i=3,4,5,6,(q=u,d,s,c,b)\hfill \\ \frac{em_b}{16\pi ^2}\overline{s}_L\sigma ^{\mu \nu }b_RF_{\mu \nu },\hfill & i=7,\hfill \\ \frac{gm_b}{16\pi ^2}\overline{s}_L\sigma ^{\mu \nu }T^ab_RG_{\mu \nu }^a,\hfill & i=8.\hfill \end{array}$$
(2.3)
Here, $`\mathrm{\Gamma }_i`$ and $`\mathrm{\Gamma }_i^{}`$ denote various combinations of the colour and Dirac matrices (see, e.g., ).

The dots in eq. (2.2) stand for UV-counterterms and non-physical operators that vanish by the QCD$`\times `$QED equations of motion. In the present section, we neglect everything that is not important for $`bs\gamma `$ at the leading order in $`\alpha _{\mathrm{em}}`$, $`m_b/M_W`$, $`m_s/m_b`$ and $`V_{ub}/V_{cb}`$. This includes other operators $`Q_i`$ of dimension 5 and 6, higher-dimensional operators, as well as terms involving leptons.

Let us assume that the decoupling of heavy particles is performed in the $`\overline{\mathrm{MS}}`$ scheme, at the renormalization scale $`\mu _0M_W`$. The values of $`C_i(\mu _0)`$ are found from the so-called matching conditions, i.e. by imposing equality of the effective- and underlying-theory Green functions at external momenta that are much smaller than masses of the decoupled particles. Next, the Wilson coefficients are evolved from $`\mu =\mu _0`$ down to $`\mu =\mu _bm_b`$, according to the Renormalization Group Equations (RGE)

$$\mu \frac{d}{d\mu }C_i(\mu )=C_j(\mu )\gamma _{ji}(\mu ),$$
(2.4)
where the anomalous dimension matrix $`\widehat{\gamma }`$ is found from UV divergences in the effective theory. This procedure results in expressing the effective Lagrangian (2.2) in terms of

$$C_i(\mu _b)=C_i^{(0)}(\mu _b)+\frac{\alpha _s(\mu _b)}{4\pi }C_i^{(1)}(\mu _b)+\left(\frac{\alpha _s(\mu _b)}{4\pi }\right)^2C_i^{(2)}(\mu _b)+\mathrm{},$$
(2.5)
where $`C_i^{(n)}(\mu _b)`$ depend on $`\alpha _s`$ only via the ratio $`\eta \alpha _s(\mu _0)/\alpha _s(\mu _b)`$. Consequently, working at a fixed order in $`\alpha _s`$, one truncates an expansion in powers of $`\alpha _s(\mu _b)`$ rather than in powers of $`\alpha _s(M_W)\mathrm{ln}(M_W^2/m_b^2)`$, as it would be the case without introduction of the effective theory. Thus, the behaviour of the perturbation series improves. This is the essence of the renormalization-group improvement in the considered case.

In the LO calculations, everything but $`C_i^{(0)}(\mu _b)`$ is neglected in eq. (2.5). At the NLO, one takes into account all the $`𝒪(\alpha _s(\mu _b))`$ contributions to $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$, including those containing $`C_i^{(1)}(\mu _b)`$.

The Wilson coefficients encode information on the short-distance QCD effects due to hard gluon exchanges between the quark lines of the leading one-loop electroweak diagrams (fig. 1). Such effects enhance the branching ratio $`\mathrm{BR}[\overline{B}X_s\gamma ]`$ by roughly a factor of three, as first pointed out in refs. .

A peculiar feature of the renormalization group analysis in $`bs\gamma `$ is that the mixing under infinite renormalization between the four-fermion operators $`Q_1,\mathrm{},Q_6`$ and the “magnetic penguin” operators $`Q_7,Q_8`$, which govern this decay, vanishes at the one-loop level. Consequently, in order to calculate the coefficients $`C_7(\mu _b)`$ and $`C_8(\mu _b)`$ at LO, two-loop calculations are necessary. Such calculations were completed in ref. . Earlier analyses  contained additional approximations or were not fully correct. The results of ref.  were subsequently confirmed in refs. .

As pointed out in refs. , the LO expression for $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$ suffers from large ($`\pm 25\%`$) renormalization scale uncertainties. Therefore, matching the experimental accuracy of eq. (1.1) requires performing a complete NLO QCD calculation. This goal has been achieved in a joint effort of many groups:

* Two-loop $`𝒪(\alpha _s)`$ corrections to the matching conditions $`C_7(\mu _0)`$ and $`C_8(\mu _0)`$ were first calculated in ref.  and subsequently confirmed by several groups .
* Two-loop mixing and one-loop matching for the four-quark operators $`Q_1,\mathrm{},Q_6`$ were found in refs. . In ref. , these results were confirmed by recalculation in a different operator basis that is more suitable for $`bs\gamma `$ analyses.
* Two-loop mixing in the sector $`(Q_7,Q_8)`$ was calculated in ref. . These results have been recently confirmed .
* Three-loop mixing between the sectors $`(Q_1,\mathrm{},Q_6)`$ and $`(Q_7,Q_8)`$ was evaluated in ref. . It is currently being verified by another group .
* The leading-order matrix elements $`s\gamma g|Q_i|b`$ and the one-loop matrix element $`s\gamma |Q_7|b`$ were calculated in refs. . Some of them were confirmed in ref.  where certain BLM corrections were included, too.
* Two-loop calculation of the matrix element $`s\gamma |Q_{1,2}|b`$ was presented in ref. . It has been recently verified and extended to the full basis of four-quark operators . The one-loop matrix element $`s\gamma |Q_8|b`$ has been found in refs. , too.

It should be emphasized that all these ingredients enter not only the analysis of $`\overline{B}X_s\gamma `$ in the SM but are also necessary in extensions of this model. The corrections to the Wilson coefficients of the operators $`Q_7`$ and $`Q_8`$ are also relevant for $`\overline{B}X_sl^+l^{}`$.

## 3 Electroweak corrections

The study of electroweak corrections begins with searching for terms that might be enhanced by large logarithms. Czarnecki and Marciano  pointed out that large logarithms $`\mathrm{ln}(m_b^2/m_e^2)`$ are absent when $`\alpha _{\mathrm{em}}^{\mathrm{on}\mathrm{shell}}`$ is used in the overall normalization of $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$.

Another type of large logarithm that might enhance some of the electroweak corrections is $`\mathrm{ln}(m_W^2/m_b^2)`$, i.e. the same logarithm that is responsible for the huge QCD enhancement of the $`bs\gamma `$ amplitude. Once $`[1\alpha _s(\mu _0)/\alpha _s(\mu _b)]0.4`$ is treated as a quantity of order unity, the considered electroweak correction is formally of order $`𝒪(\alpha _{\mathrm{em}}/\alpha _s)`$, so it might be numerically relevant, given the accuracy in eq. (1.1). However, as demonstrated in refs.  through explicit calculations, it turns out to be negligible ($`0.7\%`$).

The articles  contain results for the complete electroweak corrections to the matching conditions $`C_i(\mu _0)`$. Some of them are proportional to $`\alpha _{\mathrm{em}}(M_Z)/\mathrm{sin}^2\theta _W0.034`$. Their effect on $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$ amounts222 This number includes QED corrections to the matrix elements of $`Q_{1,2,7}`$, too. to $`1.5\%`$ for $`M_{\mathrm{Higgs}}=115`$GeV, and diminishes with increasing $`M_{\mathrm{Higgs}}`$. The authors of ref.  resolved the numerical discrepancy between refs.  and  in favor of the latter.

The only electroweak $`𝒪(\alpha _{\mathrm{em}})`$ corrections that remain unknown at present are enhanced neither by large logarithms nor by $`1/\mathrm{sin}^2\theta _W`$. Thus, we can be practically certain about their irrelevance.

## 4 Non-perturbative effects

The LO contribution to $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$ is given by the tree-level matrix element of the $`Q_7`$ operator.333 In dimensional regularization, one-loop matrix elements of $`Q_3,\mathrm{},Q_6`$ may give LO contributions, too. However, they can be absorbed into the tree level matrix element of $`Q_7`$ with a redefined Wilson coefficient . Let us temporarily assume that this operator is the only one in the effective Lagrangian (2.2), and denote the corresponding contribution to the hadronic width by $`\mathrm{\Gamma }[\overline{B}X_s\gamma ]^{(Q_7\mathrm{only})}`$.

In analogy to the analyses  of the inclusive semileptonic decay $`\overline{B}X_ue\overline{\nu }`$, one can apply the Operator Product Expansion (OPE) and Heavy Quark Effective Theory (HQET) to show that

$$\mathrm{\Gamma }[\overline{B}X_s\gamma ]^{(Q_7\mathrm{only})}=\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]^{(Q_7\mathrm{only})}\left[1+a_1\frac{\lambda _1}{m_b^2}+a_2\frac{\lambda _2}{m_b^2}+𝒪\left(\frac{\mathrm{\Lambda }_{\mathrm{QCD}}^3}{m_b^3}\right)\right].$$
(4.1)
Here, $`\lambda _{1,2}\mathrm{\Lambda }_{QCD}^2`$ are the standard HQET parameters. The value of $`\lambda _20.12\mathrm{GeV}^2`$ is known from the measured $`B`$$`B^{}`$ mass difference. The value of $`\lambda _1=(0.27\pm 0.10\pm 0.04)\mathrm{GeV}^2`$ has been determined in ref.  from the observed semileptonic $`B`$-decay spectra (see ref.  for more recent determinations). The coefficients $`a_1`$ and $`a_2`$ can be calculated within perturbation theory,444 The same refers to similar coefficients at higher orders in the ($`\mathrm{\Lambda }_{\mathrm{QCD}}/m_b`$)-expansion. which yields

$$a_1=\frac{1}{2}+𝒪(\alpha _s(m_b))\mathrm{and}a_2=\frac{9}{2}+𝒪(\alpha _s(m_b)).$$
(4.2)
The resulting $`𝒪(\mathrm{\Lambda }_{\mathrm{QCD}}^2/m_b^2)`$ non-perturbative correction on the r.h.s. of eq. (4.1) amounts to around $`3\%`$.

The relation (4.1) still holds when a lower cutoff $`E_0`$ is imposed on the photon energy in the $`\overline{B}`$-meson rest frame, provided $`E_0`$ is not too close to the endpoint $`E^{\mathrm{max}}=\frac{m_B^2m_K^{}^2}{2m_B}2.6\mathrm{GeV}`$. Acceptable values of $`E_0`$ must correspond to much larger than $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ invariant masses of the recoiling hadronic state $`X_s`$. Fig. 3 in ref.  suggests that $`E_0=1.6\mathrm{GeV}`$ is sufficiently low. More than 95% of the total $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$ originates from a peak that lays above such a cutoff.555 Consequently, $`E_0`$-dependence in eq. (4.2) can be safely neglected. This peak is now clearly seen in the $`\overline{B}X_s\gamma `$ spectrum observed by CLEO (fig. 2). Its position corresponds to the photon energy in the leading two-body decay $`bs\gamma `$.

There is neither experimental nor theoretical need to consider photons below 1.6 GeV. They are practically unobservable at the inclusive level, because of the overwhelming $`bc`$ background. On the theoretical side, keeping not too small $`E_0`$ facilitates the discussion of non-perturbative effects due to operators other than $`Q_7`$. Of course, we have to admit that 1.6 GeV is chosen arbitrarily. It could almost equivalently be 1.5 or 1.7 GeV. However, going up to the current CLEO cutoff of 2.0 GeV would increase uncertainties on the theoretical side. Data-driven extrapolation from the experimental cutoff to the theoretically preferred one is the right choice to make at present.

The discussion of non-perturbative effects becomes much more complex when we take into account operators other than $`Q_7`$. It is no longer possible to apply OPE in analogy to $`\overline{B}X_ue\overline{\nu }`$, because the $`b`$-quark annihilation and the photon emission may now be separated in space-time by more than $`\mathrm{\Lambda }_{\mathrm{QCD}}^1`$.

The contribution of $`Q_8`$ to $`\mathrm{\Gamma }[\overline{B}X_s\gamma ]`$ has been analyzed in ref.  with the help of fragmentation functions. Important non-perturbative effects have been found for low $`E_\gamma `$ only, i.e. much below $`E_0=1.6`$GeV. Thus, with our cutoff, a reliable approximation is given by the perturbative contribution to $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$ from the matrix elements of $`Q_8`$. The accuracy of this approximation does not need to be known precisely, because the perturbative contribution of $`Q_8`$ is smaller than 3%.

Similar conclusions can be drawn for the operators $`(\overline{s}\mathrm{\Gamma }b)(\overline{q}\mathrm{\Gamma }^{}q)`$, where $`q=u,d,s`$. They are present inside $`Q_3,\mathrm{},Q_6`$. Perturbative effects of their matrix elements are even smaller than that of $`Q_8`$. As far as non-perturbative effects are concerned, one might worry about production of virtual vector mesons that convert to a real photon. However, creation of such transverse mesons is impossible in the factorization approximation because $`Q_3,\mathrm{},Q_6`$ contain no $`\overline{q}\sigma _{\mu \nu }q`$ currents. Deviations from the factorization approximation are suppressed either by $`\alpha _s(m_b)`$ or by $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_b`$ . This is sufficient to make them negligible here, given the smallness of $`|C_{3,\mathrm{},6}(\mu _b)|<0.07`$, as compared to $`|C_{1,2,7,8}(\mu _b)|(0.5,\mathrm{\hspace{0.33em}1},\mathrm{\hspace{0.33em}0.3},\mathrm{\hspace{0.33em}0.15})`$.

The operators $`Q_3,\mathrm{},Q_6`$ contain $`(\overline{s}\mathrm{\Gamma }b)(\overline{b}\mathrm{\Gamma }^{}b)`$ terms, too. The $`b`$-quark loops are localized at distances much smaller $`\mathrm{\Lambda }_{\mathrm{QCD}}^1`$ in space-time. Thus, they can undergo the same treatment as $`Q_7`$, as far as non-perturbative effects are concerned. Since their perturbative contributions are minor, the non-perturbative ones are totally negligible.

Charm quark loops are the most difficult to analyze. Factorization is not sufficient here because $`2m_c/m_b`$ is not a small number. Moreover, non-factorizable contributions may be numerically important because the Wilson coefficients $`C_1`$ and $`C_2`$ are not small at all.

Let us begin with tracing down possible contributions from intermediate real $`c\overline{c}`$ states. Our cutoff $`E_0=1.6`$GeV implies that the invariant mass of the final $`X_s`$ state is smaller than $`m_{\eta _c}+m_K`$. Consequently, real $`c\overline{c}`$ states might occur only before the photon emission, i.e. in a cascade decay: $`\overline{B}Y_{c\overline{c}}X_s^{(1)}`$ followed by $`Y_{c\overline{c}}X^{(2)}\gamma `$.

The importance of such processes can be tested in the case $`Y_{c\overline{c}}=\psi `$, because separate experimental data on both (inclusive) components of the cascade decay are available. For low $`E_0`$, the resulting branching ratio of the intermediate $`\psi `$ contribution is larger than the one in eq. (1.1). It gets reduced to (a few)$`\times 10^5`$ for $`E_0=1.6`$GeV, and becomes negligible for $`E_0=2.0`$GeV .

The models used by CLEO  to extrapolate from 2.0 GeV to lower cutoffs do not include the intermediate $`\psi `$ contribution. They are based on perturbative calculations, in which the only diagram (fig. 3) that might correspond to this contribution affects $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$ by 1.7% only (for $`E_0=1.6`$GeV). Consequently, the procedure applied by CLEO is consistent with treating the intermediate $`\psi `$ contribution as background.

Identical arguments work for $`\psi ^{}`$. Higher $`c\overline{c}`$ states might produce higher energy photons. However, radiative charm annihilation processes in all the $`c\overline{c}`$ states except $`\psi `$ and $`\psi ^{}`$ have negligible branching ratios. Thus, it does not really matter whether we consider their contributions as background or not. Whatever decision is made, its effect is expected to be less than the 1.7% perturbative contribution from the diagram in fig. 3.

Having discussed the real intermediate $`c\overline{c}`$ states, we proceed to the virtual ones. Neither infrared nor collinear singularities occur in the perturbative contributions of $`c\overline{c}`$ loops to $`\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma ]`$ at NLO. Thus, according to the common wisdom, one expects that these perturbative results give reasonable estimates to the corresponding contributions to $`\mathrm{\Gamma }[\overline{B}X_s\gamma ]`$, up to corrections of order $`𝒪(\mathrm{\Lambda }_{\mathrm{QCD}}/m_{c,b})`$.

The actual situation is somewhat more complicated, because the leading one-loop diagram (fig. 4) vanishes for the on-shell photon. However, it becomes non-vanishing when a soft gluon is attached to the $`c`$-quark loop. Such a gluon may originate from the decaying $`\overline{B}`$ meson. Thus, one finds a non-perturbative effect  that is not approximated in any sense by the corresponding perturbative null. Fortunately, it can be expressed within HQET in terms of a series

$$\frac{\mathrm{\Delta }\mathrm{\Gamma }[\overline{B}X_s\gamma ]}{\mathrm{\Gamma }[\overline{B}X_s\gamma ]}=\frac{\lambda _2}{m_c^2}\underset{n=0}{\overset{\mathrm{}}{}}b_n\left(\frac{m_b\mathrm{\Lambda }_{\mathrm{QCD}}}{m_c^2}\right)^n,$$
(4.3)
in which the $`n1`$ terms are likely to be negligible, because the coefficients $`b_n`$ decrease rapidly with $`n`$ . The calculable leading $`𝒪(\mathrm{\Lambda }_{\mathrm{QCD}}^2/m_c^2)`$ term enhances the decay width by around 2.5% .

The perturbative $`𝒪(\alpha _s)`$ results described in section 2 include non-vanishing two-loop diagrams with $`c\overline{c}`$ loops, e.g. the ones obtained by adding a virtual gluon to the diagram in fig. 4. The corresponding non-perturbative effects are expected to be suppressed by both $`\alpha _s(m_b)`$ and $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_{c,b}`$. Thus, at the first glance, they might seem irrelevant. However, it remains an open question whether their suppression is numerically sufficient. No quantitative estimates of such non-perturbative effects have been performed so far. We shall discuss this issue in more detail at the end of the next section.

## 5 Phenomenological discussion

In the present section, we shall discuss the two main uncertainties in the present-day SM prediction for $`\overline{B}X_s\gamma `$. The analysis of ref.  will be largely followed.

The prediction (1.2) is obtained from the formula

$`\mathrm{BR}[\overline{B}X_s\gamma (E_\gamma >E_0)]`$ (5.1)
$`=\mathrm{BR}[\overline{B}X_ce\overline{\nu }]_{\mathrm{exp}}\left({\displaystyle \frac{\mathrm{\Gamma }[\overline{B}X_ue\overline{\nu }]}{\mathrm{\Gamma }[\overline{B}X_ce\overline{\nu }]}}\right)_{\mathrm{th}}\left({\displaystyle \frac{\mathrm{\Gamma }[\overline{B}X_s\gamma (E_\gamma >E_0)]}{\mathrm{\Gamma }[\overline{B}X_ue\overline{\nu }]}}\right)_{\mathrm{th}},`$
in which the following substitutions are made

$`\left({\displaystyle \frac{\mathrm{\Gamma }[\overline{B}X_s\gamma (E_\gamma >E_0)]}{\mathrm{\Gamma }[\overline{B}X_ue\overline{\nu }]}}\right)_{\mathrm{th}}`$ $``$ $`\left({\displaystyle \frac{\mathrm{\Gamma }[bX_s^{\mathrm{parton}}\gamma (E_\gamma >E_0)]}{\mathrm{\Gamma }[bX_u^{\mathrm{parton}}e\overline{\nu }]}}\right)_{\mathrm{NLO}}+(\begin{array}{c}\text{non-perturbative}\\ \text{corrections (}\text{4.3}\text{)}\end{array}),`$ (5.4)
$`\left({\displaystyle \frac{\mathrm{\Gamma }[\overline{B}X_ue\overline{\nu }]}{\mathrm{\Gamma }[\overline{B}X_ce\overline{\nu }]}}\right)_{\mathrm{th}}`$ $``$ $`\left({\displaystyle \frac{\mathrm{\Gamma }[bX_u^{\mathrm{parton}}e\overline{\nu }]}{\mathrm{\Gamma }[bX_c^{\mathrm{parton}}e\overline{\nu }]}}\right)_{\mathrm{NNLO}}+\left(\begin{array}{c}\text{known}𝒪\left(\lambda _2/m_b^2\right)\\ \text{corrections}\end{array}\right).`$ (5.7)
Such ratios are introduced in order to minimize uncertainties in eq. (5.1) that originate from the CKM angles and the overall factors of $`m_b^5`$. The use of $`bu`$ transitions is motivated by the fact that eq. (5.7) is known at the NNLO, while convergence of the perturbation series and non-perturbative effects are more easily controlled in eq. (5.4) than in $`\mathrm{\Gamma }[\overline{B}X_s\gamma ]/\mathrm{\Gamma }[\overline{B}X_ce\overline{\nu }]`$. The $`𝒪(\mathrm{\Lambda }_{\mathrm{QCD}}^2/m_b^2)`$ terms from eq. (4.1) have canceled in the ratio (5.4) with the analogous corrections to $`\mathrm{\Gamma }[\overline{B}X_ue\overline{\nu }]`$. On the other hand, once the quark masses are expressed in terms of the hadronic ones, the ratio (5.7) depends on both $`\lambda _1`$ and $`\lambda _2`$.

The main uncertainty in the perturbative ratio on the r.h.s. of eq. (5.4) originates from the two-loop diagrams with charm quarks presented in fig. 5. Such diagrams are the only source of $`m_c`$-dependence of this ratio. Since the higher-order (NNLO) QCD corrections are unknown, the renormalization scheme for $`m_c`$ remains arbitrary, at least within a certain class of “reasonable” schemes that do not artificially enhance the unknown corrections. As argued in ref. , the uncertainty in eq. (1.2) stemming from this scheme-dependence can be accounted for by setting $`m_c/m_b=m_c(\mu )^{\overline{\mathrm{MS}}}/m_b^{1S}`$ in the two-loop diagrams,666 Here, $`m_b^{1S}`$ stands for the $`b`$-quark mass in the so-called “1S-scheme” . It is defined as half of the perturbative contribution to the $`\mathrm{{\rm Y}}`$ mass. and varying the scale $`\mu `$ between $`m_c`$ and $`m_b`$. Such a variation is the dominant source of the error in eq. (1.2).

One could remove the considered uncertainty by calculating three-loop diagrams obtainable from fig. 5 by adding one more virtual gluon. UV-divergent parts of such diagrams have been already found in the process of calculating the NLO anomalous dimensions . Evaluating the finite parts would constitute an extremely tedious task, though not totally impossible, if numerical integration was applied. Finding the remaining NNLO corrections would be relatively simpler, given that fully automatized analytical methods are now available .

However, before undertaking such an ambitious task, one should make sure that all the non-perturbative effects are really under control. The main worry are the doubly-suppressed corrections mentioned in the last paragraph of section 4. So far, they have been neither estimated nor included in the theoretical error. They are related to precisely the same two-loop diagrams with charm quarks (fig. 5). Numerical importance of non-local parts of those diagrams777 By non-local we mean those parts that cannot be removed off-shell by finite local counterterms. can be illustrated by the fact that the r.h.s. of eq. (1.2) changes by 35% when $`m_c`$ is shifted from the original value of $`\mathrm{\hspace{0.33em}0.22}m_b`$ to the threshold for charm pair production $`m_c=\frac{1}{2}m_b`$. A $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_b`$-suppressed non-perturbative effect on the top of such a large perturbative contribution might not be negligible. Unfortunately, no systematic methods have yet been developed for calculating corrections of this type.

## 6 Conclusions

In the present paper, we have summarized the existing calculations of perturbative and non-perturbative contributions to the inclusive weak radiative $`\overline{B}`$ meson decay. We have pointed out that both the main perturbative uncertainty and the most worrisome non-perturbative effects have their origin in the fact that non-local charm quark loop contributions are particularly large. Removing the perturbative uncertainty due to $`m_c`$-dependence would be extremely tedious, but not totally impossible. However, developing a method for systematically estimating the related non-perturbative effects is desirable in advance.

The present agreement at the $`10\%`$ level between the experimental (1.1) and theoretical (1.2) determinations of $`\mathrm{BR}[\overline{B}X_s\gamma ]`$ implies that clear signatures of new physics in this observable are not likely to be found in the foreseeable future. The importance of improving the accuracy on both the experimental and theoretical sides follows from the need for strengthening the $`bs\gamma `$ constraints on beyond-SM theories. Such constraints are likely to be crucial in identifying the origin of new physics effects that we expect to encounter in the LHC era.

## 7 Acknowledgements

We would like to thank K. Baranowski, C. Bobeth, K.G. Chetyrkin, A. Czarnecki, P. Gambino, M. Jamin, A. Khodjamirian, A. Kwiatkowski, M.E. Lautenbacher, M. Münz, S. Pokorski, N. Pott, J. Urban and P.H. Weisz for fruitful collaboration in the calculations that were relevant for determining the SM prediction for $`bs\gamma `$. A.J.B. is grateful to the Max Planck Institute for Physics for financial support during the Planck02 Symposium in Kazimierz. His work presented here has been supported in part by the German Bundesministerium für Bildung und Forschung under the contract 05HT1WOA3 and the DFG Project Bu. 706/1-1. M.M. was supported in part by the Polish Committee for Scientific Research under the grant 2 P03B 121 20.


text:
# Monte Carlo Renormalization Group for Entanglement Percolation

## I Introduction

The gelation process in crosslinked polymers has been studied within the context of percolation theory Stauffer1 ; Stauffer2 ; Herrmann1 ; Herrmann2 , but the rheokinetics rheokinetics  of bulk linear polymers as a function of chain length has recieved less theoretical attention. As the density of chains, and their length increases, the viscosity starts to increase well before the onset of vitrification. Cicerone  It has been proposed deGennes ; Doi  that this is due to the entanglement of the polymer chains.

We would like to pose the question of whether the percolation of entanglement clusters is in the same universality class as percolation. We define entanglement clusters starting from ordinary sets of connected bonds. We will consider two such sets connected, as long as they come with a lattice constant of each other, i.e., share the end points of an empty bond. Since the length of the chains, or the average length between entanglements, could introduce a second length scale into the problem, this could potentially lead to a crossover to a different universality class than percolation.

It has been found in two dimensions ConiglioDaoud79  that the vulcanization process, which involves the crosslinking of long chains, is in the percolation universality class, and that there is a crossover between Self Avoiding Walk and percolation behaviour as a function of the fugacity of the crosslinkers. Similarly, Jan et al. Jan  find in three dimensions that the crossover from SAW to percolation exponents already occurs for any finite value of the concentration of initiators (from which the chainlike structures grow), for a mixture of monomers with functionalities $`2`$, i.e., again in the presence of crosslinkers. A Monte Carlo simulation in two dimensions reveals Bunde  that a growth model can cross over from SAW like behaviour to percolation, as a function of the cluster mass.

The percolation of clusters which are not necessarily connected but linked to each other by loops, has also been studied by very large MC computations Kantor . These authors find that the new critical point is very close to the ordinary percolation threshold on the cubic lattice, with the eigenvalue of the renormalized occupation probability being indistinguishable from that of the ordinary percolation problem.

In this paper we introduce a special Monte Carlo (MC) Renormalization Group procedure to investigate the universality class for the entanglement phase transition of a linear polymer system, in three dimensions and with no crosslinkers (or monomers with functionality greater than 2) present. The polymerization process is modeled by growing non-intersecting chains from a set of randomly chosen sites on a cubic lattice. Chains which occupy nearest neighbor sites on the lattice are considered to be part of the same connected cluster, and coarse grain to occupied bonds.

In the next section, we present the simulations and the MC renormalization group procedure. In section 3 we present an analysis of the results. We conclude with a discussion in section 4.

## II MC simulations and the Renormalization Group

In this section we describe a Monte Carlo renormalization group procedure, to deal with the percolation of entangled clusters of linear polymers. Since the long linear chains of the growing bulk polymer cannot be accomodated in small cells, we start by simulating the polymer growth on relatively large lattices. However we make a different choice than the one made by Swendsen Swendsen1 ; Swendsen2 ; Swendsen3  in the way that the MC renormalization group is introduced, as illustrated schematically in Fig.1.

The MC Renormalization Group procedure which we use consists of paving the large MC cell by small cells on which one performs a coarse graining transformation. Reynolds1  This induces an RG transformation on the bond occupation probability, averaged over many independent realizations.

In the conventional MC renormalization group procedure, the coarse graining is done in one step Swendsen1 ; Swendsen2 ; Swendsen3 . In the case of percolation, this would mean coarse graining the large MC cell to one occupied or empty bond - i.e., one asks the question whether the cell is spanned by a connected cluster or not. Adopting the conventional definition of a connected cluster, this one-step coarse graining procedure is clearly not appropriate in our case, since for non-intersecting linear chains, the MC cell will be spanned in case there is at least one chain that grows across it, and otherwise not. This does not describe the physical problem at hand.

The effect which we are trying to model has to do with linear chains that do not intersect, but nevertheless get physically close to each other in dense regions where they constrain each others’ motion, at least on short time scales, and thus form an effective three-dimensional network. We look, therefore, for a way in which we can capture this effectively three dimensional behaviour, and the answer lies in making a coarse graining with a scale factor $`b`$ smaller than the linear size $`L`$ of the large cell. We will illustrate below, that chains which are non intersecting can, under coarse graining, go to branched structures, which behave very much like ordinary percolation clusters.

### II.1 The Kinetic Growth Model for Linear Polymers

We have generated our configurations on a cubic lattice with linear dimension $`L`$ and periodic boundary conditions. We start from an empty lattice, and then place $`C`$ initiators at randomly chosen lattice sites. As the chains grow, the initiators are displaced in the growth direction and mark the “active end” of the chains. We don’t allow a chain to intersect itself or the other chains during this process. If the chain hits a dead end, it stops growing.

The growth rules are as follows. During one time step, from time $`t`$ to $`t+1`$, we consider each chain one by one, look at the active end of the chain and count the number of empty nearest neighbor sites. If this number is zero, the chain is trapped, and we proceed to the next chain. If the number is greater than zero, with growth probability $`k`$ we extend the chain to one of these empty sites, chosen at random. (This extra parameter $`k`$ gives us the possibility of probing arbitrarily small increments in the occupation probability on our finite cell.)We mark this new location as the active end of the chain, and place a bond between the old and the new location of this active end. After this process has been applied to all the chains, one gets the configuration at time $`t+1`$.

We use this set of kinetic growth rules, to grow a mass of non-intersecting linear chains in relatively large cells, consisting of cubic lattices of linear size $`L=8,16,32`$.

### II.2 Compactification and Coarse-Graining

In Fig. 2, we display the small cells with which our lattice will be paved, in order to perform the coarse graining transformation. This choice for the coarse graining rule has been used by Reynolds et al. Reynolds2  and Bernasconi Bernasconi . Under our coarse-graining rule each cell (with $`3b^3`$ bonds, where $`b`$ is the rescaling factor, here chosen to be 2) goes to at most three bonds in the coarse-grained lattice.

Before we coarse grain, we perform a compactification of our randomly connected clusters by connecting nearest neighbor occupied sites, which are not already connected by a bond.(See Fig. 3) This is motivated by the fact that in dense regions, chains which physically get in each others’ way are not chemically connected, but nevertheless contribute to the instantaneous shear modulus of the effective network.

We determine the state of the bonds in the coarse-grained lattice by considering the configurations in each small cell, as illustrated in Fig. 4. We look for a path which spans the cell in each Cartesian direction. If there is at least one such path this cell goes to an occupied bond in this direction in the coarse-grained lattice. Otherwise, it goes to an empty bond. By applying this rule to each small cell in the original lattice we get the new coarse-grained lattice.

It can be noticed in Fig. 4 that the occupied bonds spanning the small cell in the horizontal and the vertical directions are not chemically connected, and yet, (even without the bonds added for compactification) they go, under coarse graing to an occupied junction.

We now compute the fraction of occupied bonds in both the original and coarse-grained lattices, $`p(t)`$ and $`p^{}(t)`$, where we also count the bonds added in the compactification process. These are averaged over many independent runs, to find $`<p(t)>`$ and $`<p^{}(t)>`$, for each time $`t`$. Eliminating $`t`$, the coarse graining procedure leads to a renormalization group transformation on the average occupation probability, as shown in Figure 5. The fixed point of the RG transformation and the RG eigenvalue are computed from the MC data, leading to the value of the critical exponent we are interested in.

## III Finite size scaling and the fractal dimension

The correlation length,

$$\xi (pp^{})^\nu ,$$
(1)
for $`p`$ very close to the critical point $`p^{}`$ (the fixed point of the RG transformation) exceeds $`L`$. From then on, $`\xi L`$, i.e., it behaves like a constant with respect to $`(pp^{})`$. In other words, the relationship in Eq.(1) breaks down and a finite size scaling analysis is in order. Stauffer1

The mass contained in the incipient infinite cluster, $`\mathrm{\Delta }ML^{D_f}`$, where $`D_f`$ is the fractal dimension of the percolation cluster,so that $`P_{\mathrm{}}=\mathrm{\Delta }M/VL^{\beta /\nu }L^{D_fd}`$. We may also show that the mass $`M_0`$ contained in the rest of the finite clusters scales like $`M_0/VL^{(2\tau )D_f}`$, with $`(2\tau )D_f=D_fd`$. Thus, the total concentration $`p=(\mathrm{\Delta }M+M_0)/VM/V`$ of occupied bonds scales like

$$pL^{D_fd}.$$
(2)
We will make use of this fact to compute the fractal dimension of the percolation cluster, which can be directly related to the other known critical exponents, e.g. via $`D_f=(\beta +\gamma )/\nu `$Herrmann2

From (2), we see that under a rescaling transformation $`LL^{}`$=$`L/b`$, $`MM^{}b^{D_f}M`$, so that

$$p^{}=\frac{M^{}}{V^{}}b^{dD_f}p.$$
(3)
At the fixed point of the RG transformation, $`p^{}`$, with $`\lambda dp^{}/dp|_p^{}`$, one has, $`\lambda =b^{dD_f}`$, or

$$D_f=d\frac{\mathrm{ln}\lambda }{\mathrm{ln}b}.$$
(4)
In Fig. 5, we plot $`<p^{}(t)>`$ versus $`<p(t)>`$ for $`L`$=32 where the averages have been performed over $`10^4`$ independent runs, with $`C=0.02\times L^3`$ and $`k=0.01`$. Notice that $`p^{}0.02`$ is much smaller than the expected value for percolation on a cubic lattice (0.2488) Stauffer1 . We have varied $`C`$ between 1 to 9$`\%`$ of the lattice points, and a fixed point $`p^{}`$ has been found for all these different $`C`$ values, although the value of $`p^{}`$ descreases with decreasing $`C`$. This can be seen to result from an effective length scale introduced by the increasingly long chains growing from fewer and fewer initiators, i.e., an effective lattice spacing of the order of $`C^{1/3}/L`$ Jan .

We obtain $`\lambda `$ by making a linear fit to the data around this point, and we obtain the value of the critical exponent $`D_f`$ from Eq. (4). Our results are given in Table 1. Within our error bars we have not found any dependence of the RG eigenvalue on the concentration of initiators, within the range of values we have considered. The errors reported in Table 1 are those arising from taking one standard deviation of each of our data points reported in the inset of Fig. 5. The best value of $`D_f`$ for percolation clusters in three dimensions, obtained from simulation results Stanley  is $`2.53\pm 0.03`$. Given that we are performing a finite cell real space renormalization group calculation, our results are consistent with the percolation fractal dimension.

## IV Discussion

The physics of bulk linear polymers is an extremely interesting and rapidly growing field. As the average molecular weight (or chain length) grows, bulk linear polymers are known to exhibit many of the properties of ordinary gels, such as resistance to shear and the capacity to take up solvent and swell, while retaining their original shape. Yilmaz  This behaviour is present even in linear polymers like PMMA (Poly-methyl methacrylate) where inter-chain interactions Saiani ; Erukhimovich  are extremely weak, and therefore must arise from purely geometrical effects, in other words entanglements. Zhu ; Tian  In this context we use “entanglement” to mean a dense region where chains physically impede the motion of other chains, and not in the mathematically precise sense of knots Kantor ; Holroyd1 ; Holroyd2 ; Grimmet .

The transition to a regime where the presence of an entanglement network can be ascertained via dilatometric techniques Okay1  takes place before or around the onset of the “gel effect,” namely the characteristic rapid growth in the rate of polymerization Okay2 ; Okay3 . This effect arises due to the transition to a diffusion-limited regime, where end to end termination of chains is essentially suppressed, so that single chains grow rapidly Okay1 ; Okay2 ; Okay3 ; Zhu . Its onset has been conjectured to be due to the entanglement of the linear chains deGennes ; Doi .

Defining entanglement clusters as chains which will be considered connected if they happen to pass through nearest neighbor sites, we set out to study the scaling behaviour in the vicinity of the entanglement percolation threshold in three dimensions. We have demonstrated, via a modified MC Renormalization Group, that for realistic concentrations of initiatiors, entangled clusters of linear chains have the same fractal dimension as ordinary percolation, although the entanglement percolation threshold is much lower than for standard bond percolation on a cubic lattice.

Acknowledgements

We would like to thank Nihat Berker for useful comments. One of us would like to acknowledge partial support from the Turkish Academy of Sciences.


text:
# Magnetic behavior of single crystalline Ho2PdSi3

## Abstract

The magnetic behavior of single-crystal Ho2PdSi3, crystallizing in an AlB2-derived hexagonal structure, is investigated by magnetic susceptibility ($`\chi `$) and electrical resistivity ($`\rho `$) measurements along two directions. There is no dramatic anisotropy in the high temperature Curie-Weiss parameter or in the $`\rho `$ and isothermal magnetization data, though there is a noticeable anisotropy in the magnitude of $`\rho `$ between two perpendicular orientations. The degree of anisotropy is overall less prominent than in the Gd (which is an S-state ion!) and Tb analogues. A point of emphasis is that this compound undergoes long range magnetic ordering below 8 K as in the case of analogous Gd and Dy compounds. Considering this fact for these compounds with well-localised f-orbital, the spin glass freezing noted for isomorphous U compounds in the recent literature could be attributed to the role of the f-ligand hybridization, rather than just Pd-Si disorder.

PACS Nos. 75.50.-y; 75.30.Gw; 75.50.Ee; 75.90.+w

Corresponding author. E-mail address: sampath@tifr.res.in

The investigation of intermetallic compounds of the general formula, R2TX3 (R= rare-earth or U, T = transition metal and X = Si or Ge), derived from AlB2-type hexagonal structure , is of great interest, considering that the properties of many of these compounds are relevant to current problems in the field of magnetism and superconductivity \[see, for instance, references 1-12\]. We have also carried out intense studies on the single crystalline forms of the series R2PdSi3 for R= Ce, Gd, and Tb, and found novel anisotropic behavior in magnetic and transport data \[5-7\]. It was however noticed that such anisotropy is less evident in the case of Dy , and this finding is surprising considering that even the Gd (which is a S-state ion) compound exhibits strongly anisotropic properties. In view of this situation, we considered it worthwhile to extend single crystal investigations to other members of this series. Therefore, we have investigated single crystalline Ho2PdSi3, the results of which are reported in this article. Another motivation to take up this investigation is that the analogous U compounds have been found  to undergo spin glass freezing, whereas previous reports  for the above heavy rare-earth compounds reveal features due to long range magnetic ordering. It is therefore of interest to carefully investigate single crystals of additional f-localised systems, in order to explore possible role of f-hybridization to decide the spin glass ground state in U compounds.

The single crystals of Ho2PdSi3 (25 to 50 mm in length and 7mm in diameter) were grown by the floating zone method by rf inductive heating from stoichiometric polycrystalline feed rods. The orientation of the single crystals was determined by x-ray Laue backscattering. Pieces ( 2mm X 2mm X 10 mm) with two orientations \[10$`\overline{1}`$0\] and  were cut from this rod for measurements. Isothermal magnetization (M) measurements at 2 and 5 K and magnetic susceptibility ($`\chi `$) measurements in the temperature (T) interval 1.8 - 300 K were performed employing a commercial superconducting quantum interference device (Quantum Design) as well as a vibrating sample magnetometer (Oxford Instruments). The electrical resistivity ($`\rho `$) measurements were performed by a conventional four-probe method employing silver paint for electrical contacts.

The T dependence of $`\chi `$ at low temperatures taken at 100 Oe and 2 kOe for the zero-field-cooled (ZFC) condition of the crystal is shown in Fig. 1. For H= 100 Oe, we have in addition taken the data for the field-cooled (FC) state of the specimen. It is distinctly clear from this figure that, for H//, there is a well-defined peak at 8 K, establishing the onset of magnetic ordering. The transition temperature is higher than that reported (6 K) for polycrystals . It should be noted that there is no separation of ZFC-FC curves as in the case of single crystals of Gd, Tb and Dy alloys , which is a sufficient proof that this compound undergoes long-range magnetic ordering, rather than spin-glass freezing. As the T is lowered, there is an upturn in $`\chi `$ values below 4 K, which is indicative of possible spin reorientation effects apparently characteristic of this class of compounds . If $`\chi `$ data are collected at H= 2 kOe, this upturn vanishes; however, a broad shoulder appears around 4 K in the plot of $`\chi `$ versus T as though there is a change in the magnetic structure with decreasing T. Now turning to the data for H//\[10$`\overline{1}`$0\], the 8 K transition appears only as a broad shoulder in the plot of $`\chi `$ versus T; however, the transition below 4 K appears as a prominent feature, that is, a tendency to flatten for both the fields. The $`\chi `$ values are noticeably lower for this direction compared to those for H//. In short, $`\chi `$ behavior appears to be anisotropic in the magnetically ordered region.

With respect to the paramagnetic behavior (Fig. 2), the plot of inverse $`\chi `$ versus T is nearly linear, though there is a tendency for a small deviation below about 50 K presumably due to crystal field effects. The nature of the deviation is different for these orientations (see Fig. 2) in the sense that the experimental values at low temperatures in the paramagnetic state (say, in the range 10 - 30 K) are higher for H//\[10$`\overline{1}`$0\] with respect to the straight line extrapolated from the high temperature ($`>`$100 K) linear region, whereas for H// the values are lower. The value of the effective moment obtained from the high temperature linear region is very close to that expected for trivalent Ho ions (close to 10.6 $`\mu _B`$/Ho mol) and the Curie-Weiss temperature ($`\theta `$p) is practically the same for both directions (close to 2 K). Thus, the anisotropy is insignificant in the high temperature $`\chi `$ data. It is to be remarked that, the absolute value of $`\theta `$p is nearly the same as the Néel temperature (TN) for H//, as inferred from the linear region below 30 K; however the sign of $`\theta `$p is positive, implying the existence of significant ferromagnetic correlations. The low temperature $`\theta `$p however changes sign (to -3 K) for H//\[10$`\overline{1}`$0\] consistent with antiferromagnetic coupling.

The isothermal M behavior (Fig. 3) shows no hysteresis at 2 or 5 K. A close inspection of the data at 5 K suggests that there is a sudden increase in the slope in a field of about 2 kOe, as though there is a spin reorientation, for both the directions; the field at which this occurs increases to about 3 kOe at 2 K as shown in an expanded form in the inset for one direction. This finding suggests that the zero-field magnetic ordering is antiferromagnetic. The magnetic moment at high fields tends towards saturation to a constant value, which is indicative of field-induced ferromagnetism. The value of the saturation moment is somewhat lower than that expected for trivalent Ho ions (10 $`\mu _B`$/Ho ion) and it appears that higher fields are required to attain this saturation value. It is to be noted that, for the 2 K data for H//, the values at high fields are comparatively lower than those at 5 K, which may be a consequence of the change in spin orientation as the T is lowered. Needless to remark that, isothermal M in the paramagnetic state is a linear function of H (not shown in the figure).

In Fig. 4, we show the T dependence of $`\rho `$. The low T data has also been shown in an expanded form in the insets. Above TN, $`\rho `$ for both crystallographic directions of excitation current has a positive temperature coefficient typical of metals (as expected). Though the values for the two directions are different, the features are qualitatively the same, even below TN. The value of $`\rho `$, instead of exhibiting a drop at TN due to the loss of spin-disorder contribution, shows an upturn. This observation signals the formation of magnetic Brillouin-zone boundary gaps below 8 K, thereby endorsing the previous conclusion from neutron diffraction data  that the magnetic structure is of a modulated type (incommensurate with the lattice).

To conclude, we have reported magnetic behavior of single crystalline Ho2PdSi3 for two orientations. Among various findings, we would like to emphasize that the overall features in isothermal M, temperature dependence of $`\rho `$ and paramagnetic $`\chi `$ behavior are essentially isotropic, though the shapes of the $`\chi `$ versus T plots in the magnetically ordered state appear to show some anisotropy and the magnitude of $`\rho `$ is noticeably different for two perpendicular orientations. In this sense, the properties of this compound are similar to that of Dy2PdSi3 . In sharp contrast to this situation, the features in isothermal M and $`\rho `$ near the magnetic region are strongly anisotropic not only in the Tb analog but also for the Gd (which is an S-state ion) analog . Finally, with the observation of long range magnetic ordering in another heavy rare-earth compound of this structure in single crystalline form, it can be stated that the spin-glass behavior of the isostructural U compounds is determined in some fashion by f-hybridization (endorsing the conclusion in Ref. 12), rather than Pd-Si disorder alone, with the assumption that the degree of possible Pd-Si site disorder is the same for U as well as for the heavy rare-earth members of this series with long range magnetic order.

Some of the authors (H.B., W.L. and G.B) gratefully acknowledge the financial support of DFG by the Sonderforschungsbereich 463 ’Rare Earth-Transition-Metal-intermetallic compounds: Structure, Magnetism, Transport’.