text:
# Non–Mesonic Weak Decay of Λ–hypernuclei: a new determination of the Γ_𝑛/Γ_𝑝 ratio

## 1 INTRODUCTION

An old challenge of hypernuclear studies has been to secure the “elusive” theoretical explanation of the large experimental values ($`1`$) of the ratio, $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$, between the neutron– and proton–induced non–mesonic decay rates, $`\mathrm{\Gamma }_n\mathrm{\Gamma }(\mathrm{\Lambda }nnn)`$ and $`\mathrm{\Gamma }_p\mathrm{\Gamma }(\mathrm{\Lambda }pnp)`$ .

Because of its strong tensor component, the one–pion–exchange (OPE) model supplies very small $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ ratios, typically in the interval $`0.05÷0.20`$ for $`s`$– and $`p`$–shell hypernuclei. On the contrary, the OPE description can reproduce the total non–mesonic decay rates observed for these systems. Other interaction mechanisms are then expected to correct for the overestimation of $`\mathrm{\Gamma }_p`$ and the underestimation of $`\mathrm{\Gamma }_n`$ characteristic of the OPE. Those which have been studied extensively in the literature are the following ones: i) the inclusion in the $`\mathrm{\Lambda }NnN`$ transition potential of mesons heavier than the pion (also including the exchange of correlated or uncorrelated two–pions) ; ii) the inclusion of interaction terms that explicitly violate the $`\mathrm{\Delta }I=1/2`$ rule ; iii) the inclusion of the two–body induced decay mechanism  and iv) the description of the short range $`\mathrm{\Lambda }NnN`$ transition in terms of quark degrees of freedom , which automatically introduces $`\mathrm{\Delta }I=3/2`$ contributions.

Some progress in the theory of non–mesonic decay has been experienced in the last years. A few calculations  with $`\mathrm{\Lambda }NnN`$ transition potentials including heavy–meson–exchange and/or direct quark contributions obtained ratios more in agreement with data, without providing, nevertheless, an explanation of the origin of the puzzle . Very recently, the $`\mathrm{\Lambda }NnN`$ interaction has been studied within an effective field theory framework  with a weak decay model consisting of OPE, one–kaon–exchange and $`|\mathrm{\Delta }S|=1`$ four–fermion contact terms.

In the light of the experiments under way and/or planned at KEK , FINUDA  and BNL , it is important to develop different theoretical approaches and strategies for the determination of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ from data. In this contribution we discuss an evaluation of nucleon–nucleon coincidence distributions in the non–mesonic weak decay of $`{}_{\mathrm{\Lambda }}{}^{5}\mathrm{He}`$ and $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{12}`$C . This work is motivated by the fact that, in principle, correlation observables permit a cleaner extraction of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ from data than single–nucleon observables. This is due to the elimination of interference terms between $`n`$– and $`p`$–induced decays , which are unavoidable in experimental data and cannot be taken into account by the Monte Carlo methods usually employed to simulate the nucleon propagation through the residual nucleus. For a detailed discussion of this issue see Ref. .

The calculations are performed by combining a one–meson–exchange (OME) model describing one–nucleon induced weak decays in finite nuclei with an intranuclear cascade code taking into account the nucleon final state interactions. The two–nucleon induced channel is also taken into account, treating the nuclear finite size effects by means of a local density approximation scheme.

We also perform a weak interaction model independent analysis to extract an estimate for $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ using preliminary results from KEK  on two–nucleon angular and energy correlations. The resulting $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ values for $`{}_{\mathrm{\Lambda }}{}^{5}\mathrm{He}`$ turn out to be substantially smaller than those obtained from single nucleon distributions analyses  and fall within the predictions of recent theoretical studies .

The work is organized as follows. In Sec. 2 we give an outline of the models employed to describe the non–mesonic weak decay and we discuss the main features of the intranuclear cascade simulation accounting for the nucleon propagation inside the residual nucleus. A selection of our results is discussed in Sec. 3 and the conclusions are given in Sec. 4.

## 2 MODELS

### 2.1 Weak decay

The one–nucleon induced non–mesonic decay rates and the distributions of the nucleons produced in these processes are obtained with the OME model of Refs. . The OME weak transition potential contains the exchange of $`\rho `$, $`K`$, $`K^{}`$, $`\omega `$ and $`\eta `$ mesons in addition to the pion. The strong couplings and strong final state interactions acting between the weak decay nucleons are taken into account by using a scattering $`nN`$ wave function from the Lippmann–Schwinger ($`T`$–matrix) equation obtained with NSC97 (versions “a” and “f”) potentials . The corresponding decay rates for $`{}_{\mathrm{\Lambda }}{}^{5}\mathrm{He}`$ and $`{}_{\mathrm{\Lambda }}{}^{12}\mathrm{C}`$ are listed in Table 1 (OMEa and OMEf) together with the OPE predictions.

The differential and total decay rates for the two–nucleon induced process, $`\mathrm{\Lambda }npnnp`$, are calculated with the polarization propagator method in local density approximation (LDA) of Refs. . In such a calculation, the simple OPE mechanism, supplemented by strong $`\mathrm{\Lambda }N`$ and $`NN`$ short range correlations (given in terms of phenomenological Landau functions), describes the weak transition process. In the present calculation, the distributions of the nucleons emitted by two–nucleon stimulated decays and the value of $`\mathrm{\Gamma }_2`$ are properly scaled to maintain the ratio $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ unchanged: we then use $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1(\mathrm{\Gamma }_2/\mathrm{\Gamma }_1)^{\mathrm{LDA}}=0.20`$ for $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{5}`$He and $`0.25`$ for $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{12}`$C.

### 2.2 Intranuclear cascade simulation

In their way out of the nucleus, the weak decay (i.e., primary) nucleons continuously change energy, direction and charge due to collisions with other nucleons. As a consequence, secondary nucleons are also emitted.

We simulate the nucleon propagation inside the residual nucleus with the Monte Carlo code of Ref. . A random number generator determines the decay channel, $`n`$–, $`p`$– or two–nucleon induced, according to the values of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ and $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ predicted by our finite nucleus and LDA approaches. Positions, momenta and charges of the weak decay nucleons are selected by the same random number generator, according to the corresponding probability distributions given by the finite nucleus and LDA calculations.

Once they are produced, the primary nucleons move under a local potential $`V_N(R)=k_{F_N}^2(R)/2m_N`$, where $`k_{F_N}(R)=[3\pi ^2\rho _N(R)]^{1/3}`$ ($`N=n,p`$) is the local nucleon Fermi momentum corresponding to the nucleon density $`\rho _N(R)`$. The primary nucleons also collide with other nucleons of the medium according to free space nucleon–nucleon cross sections  properly corrected to take into account the Pauli blocking effect. For further details concerning the intranuclear cascade calculation see Ref. . Each Monte Carlo event ends with a certain number of nucleons which leave the nucleus along defined directions and with defined energies. One can then select the outgoing nucleons and store them in different ways, as we shall do in the discussion of Section 3.

## 3 RESULTS

The ratio $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ is defined as the ratio between the number of weak decay $`nn`$ and $`np`$ pairs, $`N_{nn}^{\mathrm{wd}}`$ and $`N_{np}^{\mathrm{wd}}`$. However, due to two–body induced decays and (especially) nucleon FSI effects, one expects the following inequality:

$$\frac{\mathrm{\Gamma }_n}{\mathrm{\Gamma }_p}\frac{N_{nn}^{\mathrm{wd}}}{N_{np}^{\mathrm{wd}}}\frac{N_{nn}}{N_{np}}R_2[\mathrm{\Delta }\theta _{12},\mathrm{\Delta }T_n,\mathrm{\Delta }T_p],$$
(1)
when the observable numbers $`N_{nn}`$ and $`N_{np}`$ are determined by employing particular intervals of variability of the pair opening angle, $`\mathrm{\Delta }\theta _{12}`$, and the nucleon kinetic energies, $`\mathrm{\Delta }T_n`$ and $`\mathrm{\Delta }T_p`$. The results discussed in Ref.  clearly show the dependence of $`N_{nn}/N_{np}`$ on $`\mathrm{\Delta }\theta _{12}`$ and $`\mathrm{\Delta }T_n`$ and $`\mathrm{\Delta }T_p`$. However, $`N_{nn}/N_{np}`$ turns out to be much less sensitive to FSI effects and variations of the energy cuts and angular restrictions than $`N_{nn}`$ and $`N_{np}`$ separately.

The numbers of nucleon pairs $`N_{NN}`$ —which we consider to be normalized *per non–mesonic weak decay*— are related to the corresponding quantities for the neutron– ($`N_{NN}^{1\mathrm{B}\mathrm{n}}`$) proton– ($`N_{NN}^{1\mathrm{B}\mathrm{p}}`$) and two–nucleon ($`N_{NN}^{2\mathrm{B}}`$) induced processes via the following equation:

$$N_{NN}=\frac{N_{NN}^{1\mathrm{B}\mathrm{n}}\mathrm{\Gamma }_n+N_{NN}^{1\mathrm{B}\mathrm{p}}\mathrm{\Gamma }_p+N_{NN}^{2\mathrm{B}}\mathrm{\Gamma }_2}{\mathrm{\Gamma }_n+\mathrm{\Gamma }_p+\mathrm{\Gamma }_2}N_{NN}^{\mathrm{\Lambda }nnn}+N_{NN}^{\mathrm{\Lambda }pnp}+N_{NN}^{\mathrm{\Lambda }npnnp}.$$

Here, $`N_{NN}^{1\mathrm{B}\mathrm{n}}N_{NN}^{\mathrm{\Lambda }nnn}(\mathrm{\Gamma }_n+\mathrm{\Gamma }_p+\mathrm{\Gamma }_2)/\mathrm{\Gamma }_n`$ is normalized per neutron–induced non–mesonic weak decay, etc. Therefore, $`N_{NN}^{1\mathrm{B}\mathrm{n}}`$, $`N_{NN}^{1\mathrm{B}\mathrm{p}}`$ and $`N_{NN}^{2\mathrm{B}}`$ ($`N_{NN}^{\mathrm{\Lambda }nnn}`$, $`N_{NN}^{\mathrm{\Lambda }pnp}`$ and $`N_{NN}^{\mathrm{\Lambda }npnnp}`$) do not depend (do depend) on the interaction model employed to describe the weak decay.

In Table 2 the ratio $`N_{nn}/N_{np}`$ predicted by the OPE, OMEa and OMEf models for $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{5}`$He and $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{12}`$C is given for two opening angle intervals and for $`T_n,T_p30`$ MeV. The results of the OMEa and OMEf models are in reasonable agreement with the preliminary KEK–E462 data for $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{5}`$He .

### 3.1 A weak interaction model independent analysis of data

We now discuss a weak interaction model independent analysis of KEK coincidence data. To this aim, we make use of the 6 quantities $`N_{nn}^{1\mathrm{B}\mathrm{n}}`$, $`N_{nn}^{1\mathrm{B}\mathrm{p}}`$, $`N_{nn}^{2\mathrm{B}}`$, $`N_{np}^{1\mathrm{B}\mathrm{n}}`$, $`N_{np}^{1\mathrm{B}\mathrm{p}}`$ and $`N_{np}^{2\mathrm{B}}`$ entering Eq. (3): they are quoted in Table 3 and, by definition, do not depend on the model used to describe the weak decay. We have thus to employ the following relation:

$$\frac{N_{nn}}{N_{np}}=\frac{\left(N_{nn}^{1\mathrm{B}\mathrm{n}}+N_{nn}^{2\mathrm{B}}{\displaystyle \frac{\mathrm{\Gamma }_2}{\mathrm{\Gamma }_1}}\right){\displaystyle \frac{\mathrm{\Gamma }_n}{\mathrm{\Gamma }_p}}+N_{nn}^{1\mathrm{B}\mathrm{p}}+N_{nn}^{2\mathrm{B}}{\displaystyle \frac{\mathrm{\Gamma }_2}{\mathrm{\Gamma }_1}}}{\left(N_{np}^{1\mathrm{B}\mathrm{n}}+N_{np}^{2\mathrm{B}}{\displaystyle \frac{\mathrm{\Gamma }_2}{\mathrm{\Gamma }_1}}\right){\displaystyle \frac{\mathrm{\Gamma }_n}{\mathrm{\Gamma }_p}}+N_{np}^{1\mathrm{B}\mathrm{p}}+N_{np}^{2\mathrm{B}}{\displaystyle \frac{\mathrm{\Gamma }_2}{\mathrm{\Gamma }_1}}},$$
(2)
using $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ and $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ as fitting parameters.

In Fig. 1 we report the dependence of $`N_{nn}/N_{np}`$ for $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{5}`$He on $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ for four different values of $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$. The figure corresponds to the experimentally interesting case with a nucleon energy threshold of $`30`$ MeV and the angular restriction cos $`\theta _{np}0.8`$. For a given value of $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$, Fig. 1 permits an immediate determination of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ by a direct comparison with data for the observable $`N_{nn}/N_{np}`$.

By using the $`{}_{\mathrm{\Lambda }}{}^{5}\mathrm{He}`$ data $`N_{nn}/N_{np}=0.44\pm 0.11`$ from KEK–E462  and neglecting the two–nucleon induced mechanism (i.e., $`\mathrm{\Gamma }_2=0`$), Eq. (2) supplies:

$$\frac{\mathrm{\Gamma }_n}{\mathrm{\Gamma }_p}\left({}_{\mathrm{\Lambda }}{}^{5}\mathrm{He}\right)=0.39\pm 0.11.$$
(3)
By employing the value $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1=0.2`$ (i.e., the one obtained with the model of Ref.  and used in the present calculation), a 34% reduction of the ratio is predicted:

$$\frac{\mathrm{\Gamma }_n}{\mathrm{\Gamma }_p}\left({}_{\mathrm{\Lambda }}{}^{5}\mathrm{He}\right)=0.26\pm 0.11.$$
(4)
These results for $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ are in agreement with the pure theoretical predictions of Refs. . On the contrary, they are rather small if compared with previous determinations  ($`0.93\pm 0.55`$)  ($`1.97\pm 0.67`$) from single–nucleon spectra analyses. Actually, all the previous experimental analyses of single–nucleon spectra , supplemented in some cases by intranuclear cascade calculations, derived $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ values in disagreement with pure theoretical predictions. In our opinion , the fact that our calculations reproduce coincidence data for values of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ as small as $`0.3÷0.4`$ could signal the existence of non–negligible interference effects between the $`n`$– and $`p`$–induced channels in those old single–nucleon data.

Forthcoming coincidence data from KEK and FINUDA could be directly compared with the results presented here and in Ref. . This will permit to achieve new determinations of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ —which will help in clarifying better the origin of the $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ puzzle— and to establish the first constraints on $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$.

## 4 CONCLUSIONS

To summarize, our weak interaction models supplemented by FSI through an intranuclear cascade simulation provide double–coincidence observables which are in reasonable agreement with preliminary KEK–E462 data for $`{}_{\mathrm{\Lambda }}{}^{}{}_{}{}^{5}`$He.

Through a weak interaction model independent analysis in which $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ and $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ are free parameters we reproduce the KEK $`{}_{\mathrm{\Lambda }}{}^{5}\mathrm{He}`$ data $`N_{nn}/N_{np}=0.44\pm 0.11`$ if $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p=0.39\pm 0.11`$ and $`\mathrm{\Gamma }_2=0`$ or $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p=0.26\pm 0.11`$ and $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1=0.2`$. Although these values of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ extracted from data agree with other recent pure theoretical evaluations (such an agreement has been achieved now for the first time), they are rather small if compared with the results of previous analyses from single–nucleon spectra. We suspect that non–negligible interference effects between the neutron– and proton–induced channels affected those single–nucleon analyses.

In conclusion, although further (theoretical and experimental) work is needed, we think that our investigation proves how the study of nucleon coincidence observables can offer a promising possibility to solve the longstanding puzzle on the $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_p`$ ratio.


format: text

text:
# I Introduction

## I Introduction

In recent years there has been an unprecedented rise of interest in the study of networks and their properties Newman ; Barabasi . These networks can be basically regarded as a collection of nodes linked by edges. The nodes represent individual entities, while the edges represent interaction between any pair of nodes in the network. Many natural systems, as well as engineered systems, can be represented by such networks. Some examples of such networks are found at large spatial scales, $`e.g`$, the internet or the power grid, while others are at single organismic levels, $`e.g`$, metabolic pathways Newman ; Barabasi ; Amaral . The structure and connectivity properties of such networks and the capacity and degree of connectivity of their nodes have important consequences for processes which occur on the network. Some examples of such processes are message transfer in communication networks, the spread of infectious diseases in social networks, the spread of computer viruses in computer networks and avalanche-spread in load-bearing networks Watts-newman ; Moore ; Kuperman ; Ziff ; Newman1 ; Janaki . Such network studies have been seen to have potential applications in different fields of science and technology, including epidemiology and ecology.

Several different classes of network topologies have attracted recent interest. These include small-world (SW) network models which have short paths between any two nodes and highly clustered connections Strogatz , and their variants Watts-newman ; Ziff , generalised families of small-world networks Kleinberg  and scale-free (SF) networks which incorporate the existence of ‘hubs’ viz. nodes with many more connections than the average nodeALB-Reka ; ALB-Reka-Jeong . There have also been attempts to characterise the topology of natural networks, using measures like average path lengths, clustering co-efficients, degree distributions and the existence of ‘network motifs’Milo , viz. patterns of inter-connections occurring in complex networks at numbers that are significantly higher than those in randomised networks.

It is known that the nature of the network topology significantly affects the manner in which spread processes occur on networks. A striking example of this was found by Kleinberg, where a decentralised algorithm was able to find very short paths (resulting in very short delivery times for messages) for a two dimensional network where long range connections followed the inverse square law. Another example where the topology of the network had a significant role to play is in the spread of computer viruses. It was seen that computer viruses can spread on the SF networks with zero transmission threshold Vespignani1 ; Vespignani2  and that networks of differing topologies have differing degrees of vulnerability to attack Reka-Jeong-Barabasi ; DSCall ; Cohen-Erez  and different extents of error tolerance. The importance of such studies from the point of view of applications, is obvious.

In this paper, we set up a communication network on a two dimensional lattice. Communication networks base on two dimensional lattices have been considered earlier in the context of search alogorithms Kleinberg , as well as in the context of traffic on lattices with hosts and routers Ohira ; Sole ; Fuks1 ; Fuks2 ; FN . The lattice we consider in this paper has two types of nodes, viz. regular nodes and hubs. Each regular node on the lattice has connections with its nearest neighbours, whereas some special nodes, which are designated as hubs, have connections to all the sites within a certain area of influence. Thus, nodes falling within the influence area of a hub can directly interact with the hub, or vice-versa. The degree distribution of this network is bimodal and has finite variance. We examine two distinct spread processes on this lattice network as functions of hub-density. The distribution of travel times for a message transmitted between a fixed pair of sites for this lattice shows fat fractal behaviour as a function of hub-density. If extra connections are now introduced between the hubs so that each hub is connected to two or three other hubs, the distribution crosses over to power-law behaviour. The end-to-end short cutting of hubs whose influence areas overlap also induces cross-over, but this is much milder in nature. We also examined the process of disease transfer on this network, a process that differs from message transfer in the fact that it is an undirected probabilistic process and not a directed transfer between a source and a target. Here again the transmission threshold plotted as a function of hub-density crossed over from one type of power-law behaviour to another on the introduction of two extra connections per hub.

In section II we set up the model of the communication network. In the next section (Section III) we discuss the problem of directed message transfer between source and target nodes on this network. In section IV we discuss the spread of information as a random process as in the spread of disease. The last section summarises our conclusions.

## II The model

As stated in the introduction, the proposed network is a two-dimensional lattice of nodes where every node is connected to its nearest-neighbours. A certain fraction of the total nodes are designated as hubs, viz. as nodes which have connections to all nodes within their area of influence, where the area of influence is defined as a square area around the hub, which accommodates $`(2k+1)^2`$ nodes, including the hub node, $`k`$ being the distance in lattice units (LU) from the hub to the outer most node in any principal direction within the square. All nodes within the influence area of a hub are termed its constituent nodes. A schematic representation of a hub node is shown in Fig. 1. These hubs are randomly chosen on the lattice maintaining a minimum distance, $`d_{min}`$, between any pair of hubs. This parameter ($`d_{min}`$) controls the extent of overlap between the influence regions of the hubs. See the two hubs, ‘d’ and ‘e’, in Fig. 1. The distance between them is unity, and their influence regions are maximally overlapped. As the number of hubs on the lattice goes up, the number of connections that a constituent node can have also goes up, as it increasingly acquires links with more hubs because of the overlapping of the influence regions of the hubs. Thus the model takes into account the fact that local clustering in a geographical neighbourhood can occur for many realistic models.

The hubs in this model connect to a fixed number of nodes, as in many realistic environments, the number of links that can attach to a given hub is limited by the number of ports to the hub. We examine the influence that the presence of hubs exerts on the connectivity distribution of the network. In the absence of the hubs each node has the same degree of connectivity and the degree distribution (where the degree of a node is the number of nodes it is connected to) is a delta function with a single peak at four. However, the introduction of hubs leads to completely different connectivity patterns for the network. Due to the presence of hubs, the connectivity of the nodes ranges from four (for the regular nodes) to $`(2k+1)^21`$ (for the hubs). The degree of the constituent nodes lies between the two extremes, being equal to $`C+4`$ where the constituent node lies in the overlapping influence area of $`C`$ hubs. The degree of the constituent nodes is thus a function of the hub density, which controls the extent to which the influence areas of different hubs overlap. Fig. 2 shows a typical distribution of the degree of connectivity of nodes of all types for a lattice of size $`100\times 100`$, $`k=3`$, $`d_{min}=1`$ and the hub density is $`4.0\%`$. As expected, there are two peaks in the distribution, one, a sharp one, at $`(2k+1)^21`$ corresponding to the degree of connectivity of the hubs, and the other around four (the degree of the regular sites) with a spread around the value. This spread comes from the constituent nodes of hubs with overlapping regions of influence. It is clear that the degree distribution does not fall into either of the two broad classes of networks – the SW and SF networks. It also differs from the connectivity patterns of the random graphs. (See ref ALB-Reka-Jeong  for the connectivity distributions for the SW and SF networks as well as those of random graphs, and warren ; rozenfeld  for scale-free models with local clusters and geographic separation.) The variance in the degree distribution of our model is finite unlike that in the case of scale-free networks where the variance is infinite. This fact can have very important consequences for information spread processes on a network. The importance of the variance in the spread of disease has already been noted in the case of scale-free networks. It is therefore interesting to investigate spreading processes in models like ours which do not have infinite variance, and to compare them with spreading processes seen in the case of other networks, e.g. scale free networks with and without local clustering in a geographic neighbourhood.

We now study two types of communication problems in this network viz. the problem of message transfer and that of the spread of infectious disease.

## III The problem of message transfer

We study the transfer of a message from an arbitrary source node to an arbitrary target node on the lattice. Each (ordinary) node transfers the message to the node nearest to it in the direction which will minimise the distance from the current message holder to the target. When any constituent node of a hub is the current message holder, then the node directly sends it to the hub. The lattice distance between the sender node and the hub is just one hop because of direct communication between them, thereby speeding up the process of message transmission. If the hub is the current message holder the message is forwarded to one of the constituent nodes within its influence area, the choice of constituent node being made by minimising the distance to the target. Thus the presence of hubs on the lattice increases the message transfer speed along the path, and the total travel time depends primarily on how many hubs fall on the path for the given influence radius, $`k`$. A typical travel path for such a lattice is shown in Fig. 1 and is also indicated by the label $`\mathrm{`}O^{}`$ in Fig. 3.

### III.1 Speed enhancement schemes

It is clear that the presence of hubs on the lattice makes the process of the message transfer faster than the situation when there are no hubs. However, there is further scope to enhance the speed of transmission. Here we discuss two speed enhancing schemes which can be practicably applicable to the communication process.
Scheme I: Whenever there is an overlap between influence areas of hubs, then the message is transferred from the first hub to second then to the third, so on. Intermediate constituent nodes do not participate in the process of message transfer towards the target. For example, if there is an overlap between the influence regions of two hubs, namely, A and B, along the path, then after receiving the message the first recipient constituent node of A directly forwards it to its hub, A, which directly sends it to B, from which the message is subsequently sent to the second recipient constituent node of B. The second recipient node of A and the first recipient node of B here do not take part in the message transfer. This scheme, in essence, introduces a single short-cut for two hops per pair of overlapped hubs in the travel path. A typical travel path for this scheme is indicated by the label $`\mathrm{`}I^{}`$ in Fig. 3. As we will see in the following section, the introduction of this scheme, leads to an increase in the travel speed and a reduction in the message delivery time. However, the actual reduction in travel time for a given value of $`k`$, crucially depends on the extent of overlap, which in turn depends on the minimum distance, $`d_{min}`$, between hub nodes, as well as on the hub density. e.g. for the travel path $`\mathrm{`}I^{}`$ of Fig. 3, only one pair of hubs overlap so that the reduction in travel time is just one unit.

The minimum distance between any pair of hubs, $`d_{min}`$, decides the range of overlap between the influence regions of the hubs in our network. For any choice of $`k`$ at a fixed hub density, a hub-distribution with $`d_{min}=1`$ on the lattice guarantees that the separation between any pairs of hubs is equal to or greater than $`1`$. It is easy to see that the areas of influence of a pair of hubs with $`d_{min}=1`$ have the maximum overlap. See Fig. 1. On the other hand, $`d_{min}=2k+1`$ for a given $`k`$ results in no overlapping hubs on the lattice. Scheme I can not be implemented for this case. Distributions of hubs with other intermediate values of $`d_{min}`$ have overlaps that lie between the two extremes.

Scheme II: In the second scheme to speed up the communication process we connect individual hubs with a few (in this paper, typically $`2`$ or $`3`$) other hubs selected at random. These connections, where nodes preferentially get connected to nodes having a similar degree of connectivity in the network, are called assortative connections MayAlun ; Newman2 . Under this scheme, when a hub becomes the current message holder, it first tries to send the message through one of its assortative linkages to another hub which, among all acqaintances of the hub, happens to be the nearest to the target. If the current message holder hub cannot utilise its assortative linkages because of unsuitable locations of the end-point hubs, the message is sent to the constituent node nearest to the target. A typical travel path between source $`\mathrm{`}S^{}`$ and target $`\mathrm{`}T^{}`$ is indicated by the label $`\mathrm{`}II^{}`$ in Fig. 3. Here the Manhattan distance $`D_{st}`$ between source and target is $`142`$ where the Manhattan distance is defined as $`D_{st}=|isit|+|jsjt|`$ and ($`is`$, $`js`$) and ($`it`$, $`jt`$) are the co-ordinates for the source and the target, respectively. However the travel path labelled $`\mathrm{`}II^{}`$ needs just $`50`$ steps to travel between the source and target. In comparison, the path labelled $`\mathrm{`}O^{}`$ needs $`95`$ steps and that labelled $`\mathrm{`}I^{}`$ needs $`94`$ steps.

The simulations are carried out as follows. Two nodes are selected as the source and the target at random from a lattice of a given size. The distance between them, denoted by $`D_{st}`$, is chosen to be the Manhattan distance. The number of steps required in delivery of the message from the source to the target are counted for 50 realisations of hubs for a given hub density. Then, the two nodes, $`i.e.`$, the source and the target nodes, are replaced by two other nodes selected at random from the lattice, keeping $`D_{st}`$ unchanged. Again the message transmission steps are counted for the same number of hub realisations. This is repeated for 1000 pairs of source and target nodes for a particular hub density. It should be noted that the order of averaging makes no difference. The value of $`d_{min}`$ and $`k`$ used throughout the paper are 1 and 3, respectively, unless otherwise specified.

### III.2 Cross-over behaviour

The average travel time between a source and target a fixed distance apart is a good measure of the efficiency of the network for message transmission. This clearly depends on the density of hubs in the network, as well as on the way in which these hubs are connected. We study the behaviour of average travel times as a function of hub density for a fixed Manhattan distance $`D_{st}`$ between source and target. Our simulations were carried out for a lattice of $`500\times 500`$ nodes $`D_{st}=712`$ for the original networks as well as the networks modified by Schemes I and II. Fig. 4 shows the dependence of the average travel times, $`t_{avg}`$, as a function of hub density for the original networks (diamonds) for $`d_{min}=1`$. The average travel times decrease exponentially as the hub density increases. The data can be fitted well by the exponential function FN1 , $`f_1(x)=Q_1\mathrm{exp}[A_1x^{\alpha _1}]`$ where $`\alpha _1=0.4482`$, $`A_1=0.0142`$ and $`Q_1=730`$. This can be rewritten in the form $`f_1(x)=Q_1\mathrm{exp}[(\frac{x}{x_0})^\alpha ]`$ where $`x_0`$ is $`(\frac{1}{A_1})^{\frac{1}{\alpha _1}}`$ and has the approximate value 13259. Expanding this we get $`f_1(x)=Q_1(1X+\frac{X^2}{2!}\mathrm{})`$ where $`X=(\frac{x}{x_0})^{\alpha _1}`$. Retaining terms to the lowest order, we see that the dependence of average travel times on hub density is given by $`t_{avg}Q_1(1a\rho _{hub}^{\alpha _1})`$, an instance of fat fractal like behaviour. We plot the average travel times as a function of hub density on the same plot for the Scheme I and Scheme II networks. The behaviour of $`t_{avg}`$ in the case of Scheme I networks (plus signs) is slightly different from that of the original network case, and the exponential function acquires a mild power-law correction. The data for Scheme I networks is fitted well by a function $`f_2(x)=Q_2\mathrm{exp}[A_2x^{\alpha _2}]x^\delta `$ where $`\alpha _2=0.46`$, $`A_2=0.0145`$, $`Q_2=735`$ and $`\delta =0.00005`$. Scheme II networks show distinctly different behaviour. The Scheme II network data (boxes for two assortative connections and crosses for three assortative connections) can be fitted by a function, $`g(x)=Sx^\beta `$ where $`\beta =0.2`$ and $`S`$ is a positive constant and thus show power-law behaviour. We note that the same power law fits both sets of Scheme II data.

Fig. 5 shows the same data as Fig. 4 on a log-log plot. Here the drastic difference seen for the behaviour of average travel times in the case of networks with Scheme II operational compared to the other two networks can be clearly seen. The log-log plot of $`t_{avg}`$ against the hub density is a straight-line with slope $`\beta =0.2`$. Thus $`t_{avg}\rho _{hub}^\beta `$ shows power-law behaviour. It is again clear that the same power-law is seen for the two Scheme II cases. Thus the addition of a very small number of assortative connections per hub has induced a cross-over to power-law behaviour from the fat fractal behaviour seen in the non-assortative cases. The rate of fall off of average travel times with increase in hub density is much faster for the assortative case even in the case of low hub densities. Thus the addition of assortative connections can increase the communication efficiency of networks without increasing the number of hubs in the network especially at low hub densities.

We note that the cross-over from fat-fractal behaviour to power-law behaviour is insensitive to perturbations of the regular lattice geometry. We verified that upto $`10\%`$ variation in the $`x`$ and $`y`$ location in the location of the hub-nodes makes no difference to the cross-over behaviour, although the numerical values of the exponents change. Similar perturbations to the location of all the nodes in the lattice also make no difference to the existence of the cross-over. Changing the shape of the influence areas from squares to circular regions introduces a variation in the degree of each hub, and also a perturbation to the over-all degree distribution. We note that the cross-over is robust to this change, as well. Fig. 5 (c) show the behaviour of the hub density for the lattice network with a perturbed node distribution with a circular influence area for the hubs, with all other parameters as in Fig. 5(a) and Fig. 5(b). It is clear that the cross-over is completely stable to perturbation. The travel times can be fitted by the functions $`F(x)=635\mathrm{exp}(0.0104x^{0.417})`$ for the perturbed lattice with no assortative connections, and by the function $`G(x)=735x^{0.165}`$ for the perturbed case with assortative connections. We note that the exponent $`\beta `$ now takes the value $`0.165`$ from the value $`0.2`$ seen earlier, and the power-law $`\alpha _1`$ changes to $`0.417`$ from the value $`0.4482`$ but the constant $`A_1`$ barely changes from $`A_1=0.0142`$ to $`A_1=0.0104`$ a change in the third place after the decimal. Thus the fat-fractal behaviour of the type $`C(1x_1^\alpha )`$ remains unaltered due to the perturbation as does the power-law behaviour for the network with assortative connections. Thus, the cross-over is robust to perturbation.

### III.3 Distribution of travel times

The probability distributions of travel times on the network for the three cases above can be seen in Figs 6 and 7. It is clear that as hub densities increase, the peak shifts towards short travel times for the original network. When the message was transmitted using Scheme I, the distribution changes very slightly with virtually no effect at low hub densities. However, when Scheme II was operational, the distribution showed a marked change (See Fig. 7). The distributions of Fig. 6 (the original lattice and Scheme I) are sharply peaked about a mean travel time and are symmetric about the mean whereas those of Fig. 7 show a much wider spread and are skewed. There is also a difference between the two distributions as shown in Fig. 7. At low hub densities the distribution is bimodal, whereas the bimodality is smoothed out at higher hub density. The right peak is reminiscent of the distribution for a lattice without any hub-to-hub connections at low hub density. Overall, the new distribution indicates that even sparse hub-to-hub connections are quite capable of inducing short paths at all hub densities. The distribution also shifts to lower values of travel times demonstrating the success of the short-cutting assortative strategy (Scheme II).

## IV Communication networks vs. spread of infections

The message transfer problem discussed above involves the consistent directed transfer of a message towards a target. Each temporary message holder transfers the message, with probability one, in the direction which takes the message towards the target. Other types of processes of information spread such as the spread of computer viruses, infectious disease, rumours, popular fashions etc. have distinctly different mechanisms of spread. These processes are not directed processes, and incorporate stochastic elements in the mechanism of spread. Hence communication networks which may be very efficient for information dissemination of the first type, may not be at all, or less suitable, for the second kind. In this section we examine the spread of an infection in a population of susceptible individuals for our network. The population contains both infected and susceptible individuals. The individuals constitute the nodes of the network, and social interactions among them constitute the edges by which infection can be transferred from one node to another. A similar model is also applicable to the spread of computer viruses. We consider a single point entry for infection in our study.

Many recent studies have focussed, on this kind of spread of diseases on networks Watts-newman ; Moore ; Kuperman ; Ziff ; Newman1 ; Vespignani1 ; Vespignani2 . A susceptible individual can get infected (with some probability) only when he directly or indirectly encounters infectious individuals in the population. The structure of the contact network has important implications for three things. The first is the rate at which an infection can spread across the network, the second is the transmission threshold, $`i.e.`$ the smallest probability of infection, with which the infection can spread to a significant fraction of the nodes of the network, and the third is the choice of an effective immunization strategy.

Recent studies of disease spreading viruses or computer virus spreads on networks of the SF type show that the possibility of these viruses being persistent in the population, and of their being resurrected causing repeated epidemics, is almost independent of any transmission threshold Vespignani1 . Thus, in the case of the SF networks, all that is needed for the infection to spread throughout the population, is the occurrence of a single-time point entry of the infection into the population. Since it has been argued that in the case of sexually transmitted diseases, particularly that of the HIV/AIDS infection, the underlying contact networks have scale-free character Liljeros , this result also has implications for the spread of such diseases. However, the existence of a vanishing threshold in the case of human disease does not fit well with the conventional epidemiologist’s view-point. The infinite variance of the node-connectivity in the SF network, has been identified as the causative agent for the vanishing threshold seen in such networks MayAlun ; AlunMay ; Victor . Thus, though traditional epidemiology acknowledges the importance of heterogeneity in the rapid spatial spread of diseases like HIV/AIDS, the extreme heterogeneity seen in SF networks may not make them good candidates for modelling other kinds of infection in human or animal populations, where social contacts generally do not conform fully with the SF characteristic Amaral ; Holme . Hence the study of disease spread on other types of networks such as ours, which have finite variance, is important.

As mentioned above, network topologies also determine the choice of effective immunization strategies. Recent studies on immunization strategies for efficient and successful controlling of epidemics on SF networks have unanimously advocated the immunization of the most connected nodes MayAlun ; Zoltan ; Vespignani3 . The choice of immunization strategy would have to be quite different in the case of networks where the node-connectivity distribution is not of the SF type. Even in the case of SF networks, immunization strategies which concentrate on immunizing the most connected nodes do not take into account the finite probability that a few infected nodes, which may not be the most connected nodes on the network, might have long range connections to new regions of the susceptible population, and hence, could transmit the disease to regions quite distant from their place of origin where a fresh epidemic can ensue. Our network, which incorporates geographic distance and local clustering, is useful for such a study.

We studied information spread via an SIR type of infection process on this network. The process has the following features: (i) all nodes are equally susceptible, (ii) infection always starts from a single site, (iii) an infectious node can infect any of the nodes it is connected to with probability $`p`$, (iv) this transmission probability is the same $`i.e.`$ it has the value $`p`$, irrespective of whether the infected node is a hub or an ordinary node, (v) none of the nodes get infected twice. Also once a node becomes infected, it remains infected until it infects any one of its neighbours. Here we assume that within a typical infectious period an infected node certainly infects, at least, one susceptible node. In this study we restrict ourselves to a single episode of epidemic break.

We study the manner in which the infection spreads on a $`100\times 100`$ node lattice. Fig. 8 shows the manner in which the infection wave-front from a single infected node travels on a lattice with no hubs at $`t=5`$, $`t=10`$ and $`t=15`$ (top row, labelled A), on a lattice with hubs (middle row labelled B) at the same time steps , and a lattice with hubs with two assortative connections (bottom row labelled C). The rapid spread of infection in the bottom row can be very clearly seen. All the results in this section study infection spread for a $`100\times 100`$ lattice with $`d_{min}=1`$ and $`k=3`$. Other parameters are given in the figure captions. Each run is taken for $`400`$ time-steps.

### IV.1 Threshold behaviour and cross-over

The threshold value of the transmission probability, $`p_{th}`$, is a very crucial factor in the spread of disease. This quantity may depend on the structure and topology of the network. In the case of the spread of computer viruses on scale-free networks, the threshold value turns out to be zero, leading to very rapid spread of virus across such networks. In the case of immunological diseases, the threshold probability is finite. In this section we examine the behaviour of this quantity for our networks as a function of the hub density.

The threshold transmission probability, $`p_{th}`$, is defined in the following way. For a fixed run of $`t`$ time steps, the threshold value is defined to be the smallest value of the transmission probability for which at least 50$`\%`$ of the total number of susceptible sites are infected by half the run. Fig. 9 shows the behaviour of $`p_{th}`$ against the number of hubs. The diamonds show the observed behaviour for the original network and the plus signs the behaviour for the Scheme II network when there are two assortative (hub-to-hub) connections per hub. The data for the original network could be fitted using a function, $`F(x)=D\mathrm{exp}[Mx^\zeta ](x+1)^{\xi _1}+C_1`$. The behaviour for the Scheme II networks (plus signs) was somewhat different and had to be fitted using another function, $`G(x)=D\mathrm{exp}[Mx^\zeta ](x+1)^{\xi _2}+C_2`$, where $`D=0.25`$, $`M=0.0054`$, $`\zeta =0.975`$, $`\xi _1=0.065`$, $`\xi _2=0.13`$, $`C_1=0.046`$ and $`C_2=0.036`$. Thus, there is a change in the scaling behaviour of the threshold probability between the original network and the scheme II case. However this cross-over is a gentle change from one power-law to another, unlike the drastic change seen for the behaviour of the average travel times for the same case. The change is more pronounced at low hub density as in the other case. The reason as to why the cross-over is drastic in the case of the message transfer process but is more gentle in this case of the infection spread, may lie in the fact that in the case of disease spread, the variance appears to be the crucial quantity which determines the threshold behaviour. We plot the variance as a function of hub-density in Fig. 10 for the original network (diamonds) as well as the network with assortative connections (plus signs). It is clear that there is very little difference between the behaviour of the variance in the two cases. This appears to be the reason for the gentle cross-over.

### IV.2 Immunization strategies

Studies of the SF networks have emphasised the importance of high connectivity nodes, and immunisation strategies which immunise nodes of high connectivity are the most successful. However, nodes of moderate connectivity can open up new regions to the spread of infection if they are connected to nodes which are far away. Real-world epidemic events have always been effected by the long-distance movements of causative agents into susceptible regions Mary . Some recent studies Ahmed ; Anderson  provide evidence for this. While studies of the SW networks take cognizance of the fact that that long range connections are important in the spread of information, it is never clear what fraction of connections in a small world are really long range, since the network is always stochastically generated using re-wiring probabilities. Local clusters also play a crucial role in the spread of infection.

We try to isolate the contribution of long range connections and that of the local clustering property in information and disease spread on the network with two assortative connections per hub by the use of different immunisation strategies. Fig. 10 plots the number of new infections as a function of time for this network. The plus signs indicate the the number of new infections for the unimmunised network. The first immunisation strategy immunises the bonds which connect hubs separated by a distance greater than or equal to $`100`$ lattice units (this is the Manhattan distance between the two hubs). For this lattice the average number of bonds of this type is about $`17\%`$ of the total number of bonds. However, the hub itself is not immunised and can infect the local cluster. It is clear that the immunisation of the long range bonds causes the number of new infections to decrease (see plot with diamonds), and also cause the infection to spread more slowly. However, this effect is not pronounced. The second immunisation strategy innoculates both the hubs which are connected by such long-range bonds so that no infection travel from the hub to any of the nodes connected to it. The number of new infections for this case is plotted with boxes in the same figure. It is clear that the rate of spread of infection slows down, the number of new infections peaks at a much lower value, and the distribution develops a long tail. The last plot, viz. the plot with crosses shows the number of new infections as a function of time if the same number of hubs is immunized as in the last strategy, but if the hubs are randomly chosen (i.e. the hubs chosen do not necessarily have long range bonds). It is interesting to note that there is very little difference between the two distributions. Thus, the local cluster which connects to other local clusters (not necessarily distant ones) appears to play a more crucial role in the spread of infection than the existence of long range bonds. This is unlike the behaviour seen earlier in the case of small world networks. Thus immunisation strategies which target arbitrary local clusters are as successful as those which target local clusters with additional long range bonds. We observed this behaviour at both low and high hub densities.

## V Conclusions

To summarise, in this paper we have studied information spread on a two-dimensional communication network with nodes of two types, ordinary nodes which are connected to their nearest neighbours and hubs which are connected to all nodes within a certain range of influence. The degree distribution for this lattice is bimodal in nature, and has finite variance. The average travel time for directed message transfer between source and target on this lattice shows fat fractal behaviour as a function of the hub density, however the introduction of a small number of assortative connections between the hubs (Scheme II) induces a cross-over to power-law behaviour for this average travel time. In the case of Scheme I networks, where a short-cut was introduced between end to end hubs for consecutively overlapped hubs, a much milder cross-over was seen. We also study the spread of infection on this network by the SIR process. The threshold level for the infection probability is finite for the networks with and without assortative connections, due to the fact that both networks have finite variance. However, the threshold level as a function of hub density shows cross-over behaviour when assortative connections are introduced when compared with the original network. However, this cross-over is gentle in comparison to that observed for the average travel times for the directed message transfer for the same Scheme II case. Thus, while network topology modifies the way in which information spreads on a network, the effect appears to be stronger for directed processes than for undirected processes. We also study the spread of infection and immunisation strategies for this network, and conclude that local clustering plays as important a role as the existence of assortative connections in the rate of spread of infection. Thus assortative connections play a more crucial role in message transfer processes than in the spread of infection.

Our results can be of practical utility in a variety of contexts. In the case of directed message transfer, at low values of hub density, the average travel time between source and target can be reduced very rapidly by the introduction of very few assortative connections per hub. This is a very efficient way of reducing travel time without the introduction of new hubs. Long range connections between hubs cut the travel time drastically in these cases. On the other hand, the existence of a local clusters which can connect to other local clusters (not necessarily distant ones) seems to play an important role in the spread of infection. Thus immunization strategies which target local clusters appear to be called for. It is thus important to note that different elements of the network topology appear to be important for different types of information spread processes. We hope to explore this direction further in future work.

## VI Acknowlegement

We thanks CSIR, India for partial support for this work.


text:
# 1 Introduction

## 1 Introduction

It is becoming increasingly clear that unstructured proteins play an important biological role . In many cases, such proteins adopt a specific structure upon binding to their biological targets. Recently, it was demonstrated that the in vitro evolved ZSPA-1 protein  exhibits coupled folding and binding .

ZSPA-1 is derived from the Z domain of staphylococcal protein A, a 58-amino acid, well characterized  three-helix-bundle protein. ZSPA-1 was engineered  by randomizing 13 amino acid positions and selecting for binding to the Z domain itself. Subsequently, the structure of the Z:ZSPA-1 complex was determined both in solution  and by crystallography . In the complex, both ZSPA-1 and the Z domain adopt structures similar to the solution structure of the Z domain. However, in solution, ZSPA-1 does not behave as the Z domain; Wahlberg et al.  found that uncomplexed ZSPA-1 lacks a well-defined structure, and that its melting behavior is less cooperative than that of the Z domain.

The Z domain is a close analog of the B domain of protein A, a chain that is known to show two-state folding without any meta-stable intermediate state . The folding behavior of the B domain has been studied theoretically by many different groups, including ourselves, using both all-atom  and reduced  models. In many cases, it was possible to fold this chain, but to achieve that most models rely on the so-called Gō prescription . Our model  is, by contrast, sequence-based. This makes it possible for us to study both ZSPA-1 and the wild-type Z domain and compare their behaviors, using one and the same model.

The purpose of this note is twofold. First, we check whether our model can explain the difference in melting behavior between ZSPA-1 and the wild-type sequence. Second, using this model, we study the structure of ZSPA-1.

## 2 Materials and Methods

### 2.1 Model

The model we study  is an extension of a model with three amino acids  to a five-letter alphabet. The five amino acid types are hydrophobic (Hyd), polar (Pol), Ala, Pro and Gly. Hyd, Pol and Ala share the same geometric representation but differ in hydrophobicity. Pro and Gly have their own geometric representations.

The Hyd, Pol and Ala representation contains six atoms. The three backbone atoms N, Cα and C and the H and O atoms of the peptide unit are all included. The H and O atoms are used to define hydrogen bonds. The sixth atom is a large Cβ that represents the side chain. Gly lacks the Cβ atom but is otherwise the same. The representation of Pro differs from that of Hyd, Pol and Ala in that the H atom is replaced by a side-chain atom, Cδ, and that the Ramachandran angle $`\varphi `$ is held fixed at $`65^{}`$.

The degrees of freedom of our model are the Ramachandran torsion angles $`\varphi `$ and $`\psi `$, with the exception that $`\varphi `$ is held fixed for Pro. All bond lengths, bond angles and peptide torsion angles ($`180^{}`$) are held fixed.

The interaction potential

$$E=E_{\text{loc}}+E_{\text{ev}}+E_{\text{hb}}+E_{\text{hp}}$$
(1)
is composed of four terms. The first term is a local $`\varphi ,\psi `$ potential. The other three terms represent excluded volume, backbone hydrogen bonds and effective hydrophobic attraction, respectively (no explicit water). For simplicity, the hydrophobicity potential is taken to be pairwise additive. Only Hyd-Hyd and Hyd-Ala Cβ pairs experience this type of interaction. In particular, this means that Ala is intermediate in hydrophobicity between Hyd and Pol. The amino acids in the Hyd class are Val, Leu, Ile, Phe, Trp and Met, whereas those in the Pol class are Arg, Asn, Asp, Cys, Gln, Glu, His, Lys, Ser, Thr and Tyr. A complete description of the model, including numerical values of all the parameters, can be found in our earlier study .

In this earlier study, the model was applied to the 10–55-amino acid fragment of the B domain of protein A. Despite the simplicity of the potential, this sequence was found to have the following properties  in the model:

* It does make a three-helix bundle with the native three-helix-bundle topology,There are two possible three-helix-bundle topologies; if we let the first two helices form a U, then the third helix can be either in front of or behind this U. although the suppression of the wrong topology is not very strong. All helices are right-handed, as they should.
* Energy minimization restricted to the thermodynamically favored (native) topology gives a structure with a root-mean-square deviation (RMSD) of 1.8 Å from the NMR structure  (calculated over all backbone atoms).
* The collapse transition is much more cooperative for this sequence than for random sequences with the same composition. Moreover, chain collapse and helix formation occur at approximately the same temperature.

The relative order of chain collapse and helix formation depends strongly on the relative strength of the hydrogen bonds and the hydrophobic attraction, so the last conclusion may seem somewhat arbitrary. However, the chain does not fold to a helical bundle if the hydrogen bonds are too strong, and it does not fold in a cooperative manner if the hydrogen bonds are too weak . As a result, with our ansatz for the potential, there is not much freedom left in the choice of these parameters, if the chain is to fold to a compact helical bundle in a cooperative manner.

In the present study, we apply the same model, with unchanged parameters, to ZSPA-1 and the Z domain of protein A. Following previous calculations for the B domain , we consider the 9–54-amino acid fragments of these two sequences (corresponding to the 10–55-amino acid fragment of the B domain). It should be mentioned that we also performed calculations for the 4–54-amino acid fragments of ZSPA-1 and the Z domain, with similar results.

The amino acid sequences of ZSPA-1 and the Z domain differ at 13 positions, all of which are located in the section 9–35. Table I shows this part of the sequences.

### 2.2 Numerical Methods

To simulate the thermodynamic behavior of this model, we use simulated tempering , in which the temperature is a dynamic variable. Details on our implementation of this method can be found elsewhere . For a review of simulated tempering and other generalized-ensemble techniques, see .

In conformation space we use two different elementary moves: first, the pivot move in which a single torsion angle is turned; and second, a semi-local method  that works with seven or eight adjacent torsion angles, which are turned in a coordinated manner. The non-local pivot move is included in our calculations in order to accelerate the evolution of the system at high temperature, whereas the semi-local method improves the performance at low temperature.

Our simulations are started from random configurations. All statistical errors quoted are 1$`\sigma `$ errors obtained by analyzing data from eight independent runs.

The temperatures studied range from $`0.87T_\text{m}`$ to $`1.43T_\text{m}`$, $`T_\text{m}`$ being the melting temperature for the wild-type Z domain. The experimental value of this temperature is $`T_\text{m}=75^{}`$. Hence, the lowest and highest temperatures studied correspond to $`31^{}`$C and $`225^{}`$C, respectively. In the dimensionless energy unit used in our earlier study , $`T_\text{m}`$ is given by $`kT_\text{m}=0.630\pm 0.001`$, $`k`$ being Boltzmann’s constant. In the model we define $`T_\text{m}`$ as the maximum of the specific heat.

## 3 Results and Discussion

Using the model and methods described in the previous section, we study the 9–54-amino acid fragments of ZSPA-1 and the wild-type Z domain. The latter sequence differs only by a one-point mutation from the previously studied 10–55-amino acid fragment of the B domain. Our results for the Z domain are similar to those for the B domain  summarized in the previous section. Figure 1 shows the free energies $`F(\mathrm{\Delta },E)`$ and $`F(\mathrm{\Delta },E_{\text{hb}})`$ for the Z domain at $`T=0.87T_\text{m}`$, where $`\mathrm{\Delta }`$ denotes RMSD from the NMR structure  (PDB code 2SPZ, model 1). Two major minima can be seen, with similar hydrogen-bond energies. These minima correspond to the two possible three-helix-bundle topologies. Both topologies are significantly populated, but the average total energy is slightly lower for the native topology, and this topology is the thermodynamically favored one. We also performed an energy minimization for the native topology, by applying simulated annealing combined with a conjugate gradient method to a large number of low-temperature conformations. The minimum-energy structure obtained this way is schematically illustrated in Fig. 2. It has an RMSD of $`\mathrm{\Delta }=1.7`$ Å from the NMR structure. The corresponding result for the B domain was, as mentioned earlier, $`\mathrm{\Delta }=1.8`$ Å.

Let us now compare the behavior of the Z domain with that of the engineered ZSPA-1 sequence. By CD, Wahlberg et al.  found ZSPA-1 to be less helical than the wild-type Z domain, the mean residue ellipticity for ZSPA-1 being 60% of that for the wild-type sequence. Furthermore, they found that the helix formation sets in at a lower temperature and is less cooperative for the engineered sequence. Figure 3a shows the helix content against temperature in our model, for both sequences.We define helix content in the following way. Each amino acid, except the two at the ends, is labeled h if $`90^{}<\varphi <30^{}`$ and $`77^{}<\psi <17^{}`$, and c otherwise. The two amino acids at the ends are labeled c. An amino acid is said to be helical if both the amino acid itself and its nearest neighbors are labeled h. The total number of helical amino acids is denoted by $`N_\text{h}`$. The maximum value of $`N_\text{h}`$ is $`N4`$ for a chain with $`N`$ amino acids. In agreement with the experimental results, we find that ZSPA-1 has a lower helix content, and that the helix formation is shifted toward lower temperature for this sequence. Figure 3b shows the temperature dependence of the radius of gyration. We find that ZSPA-1 is more compact than the Z domain. A comparison with Fig. 3a shows that, in our model, chain collapse occurs before helix formation for ZSPA-1. From Fig. 3 it can also be seen that the melting behavior is less cooperative for ZSPA-1 than for the Z domain. This conclusion is supported by our data for the specific heat (not shown); the peak in the specific heat is more pronounced for the Z domain than for ZSPA-1.

That the model predicts ZSPA-1 to be more compact than the Z domain is not surprising, given that the number of hydrophobic amino acids is larger for ZSPA-1 (14) than for the Z domain (11). In addition, ZSPA-1 has one more Pro than the wild-type sequence, which does change the local properties of the chain and could affect the overall size, too. It should be pointed out that the effect of a Pro on the overall size may be poorly described by the model because all peptide bonds, also those preceding a Pro, are held fixed (trans).

The reduced helix content of ZSPA-1 shows that this sequence does not make a perfect three-helix bundle, but does not tell how its structure differs from a three-helix bundle. It could be that one of the three helices is missing and that the other two are still there, but it could also be that the disorder is more uniform along the chain, so that all three helices are present but partially disordered. The NMR analysis of ZSPA-1  does not exclude any of these two possibilities.

Figure 4 shows how the helix content varies along the chains in our model. A comparison with experimental data for the Z domain  shows that the first half of helix II is somewhat distorted in the model. As a result, it is possible that the model underestimates the structural change produced by the mutation Glu25Pro (see Table 1), which should have a helix-breaking effect. Our results for helices I and III of the Z domain are, by contrast, in good agreement with experimental data. These two helices respond very differently to the mutations leading to ZSPA-1. Our results suggest that helix III, which itself is free from mutations, remains stable in ZSPA-1, whereas helix I, which contains seven mutations (see Table 1), turns unstable. Two possible explanations why the model predicts helix I to become unstable are that the hydrophobicity pattern of helix I is less helical in ZSPA-1, and that one of the mutations, Phe13Gly, increases the flexibility of this part of the chain.

To further investigate how the mutations affect different parts of the chain, we also perform an RMSD-based analysis. For each conformation, we compute two RMSD values, $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$, for the first and second halves of the chain, respectively. The two parts of the chain are separately superimposed on the NMR structure. Figure 5 shows the probability distributions of $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$ both for ZSPA-1 and the Z domain. In line with the results in Fig. 4, we find that the two $`\mathrm{\Delta }_2`$ distributions are similar, although the distribution for ZSPA-1 is slightly wider. The two $`\mathrm{\Delta }_1`$ distributions differ, by contrast, markedly, the mean being significantly higher for ZSPA-1 than for the wild-type sequence. These results confirm that, in our model, the disorder of ZSPA-1 is not uniformly distributed along the chain; the main difference between ZSPA-1 and the Z domain lies in the behavior of the first half of the chain.

## 4 Conclusion

Using a model that combines a relatively detailed chain representation with a simple interaction potential, we have studied the thermodynamic behaviors of an engineered sequence and its parent. The model is sequence-based, which makes it possible to compare the two sequences in a straightforward manner. Despite the simplicity of the potential, we found that the model is able to capture important effects of the mutations; the mutated sequence, ZSPA-1, shows a reduced helix content and a melting behavior that is less cooperative, compared with the wild-type sequence. We also found that chain collapse occurs before helix formation sets in for ZSPA-1, and that the main difference between the two sequences lies in the behavior of the first half of the chain, which is less stable for ZSPA-1. To decide whether or not these two predictions are correct requires further experimental data.

Acknowledgments

We thank Torleif Härd for a helpful discussion. This work was in part supported by the Swedish Research Council.


text:
# Slowing heavy, ground-state molecules using an alternating gradient decelerator

## Abstract

Cold supersonic beams of molecules can be slowed down using a switched sequence of electrostatic field gradients. The energy to be removed is proportional to the mass of the molecules. Here we report deceleration of 174YbF, which is 7 times heavier than any molecule previously decelerated. We use an alternating gradient structure to decelerate and focus the molecules in their ground state. We show that the decelerator exhibits the axial and transverse stability required to bring these molecules to rest. Our work significantly extends the range of molecules amenable to this powerful method of cooling and trapping.

It has recently become possible to decelerate cold supersonic beams of light polar molecules in suitable excited states using electrostatic forces. This has been demonstrated with beams of CO Bethlem(1)99 , ND3 Bethlem(1)00 ; Bethlem(1)02 , and OH Bochinski(1)03 . Once the molecules are sufficiently slow they can be trapped electrostatically Bethlem(1)02 ; Crompvoets(1)01 . The key feature of an electrostatic decelerator is that molecules move through the fringe field of a capacitor from low to high electric field. If the molecules are weak-field-seekers the increase in their potential energy causes them to slow down. The field is then switched off and the molecules fly out of the capacitor without being accelerated. With many capacitors along the beam line, this process can be repeated until the molecules are brought to rest. The highest electric fields are on the surfaces of the electrodes, resulting in a transverse force towards the beam axis that focuses the weak-field-seeking molecules as they propagate through the decelerator Bethlem(1)02 .

Until now, no ground state beam has been decelerated because ground state molecules are strong-field seekers. Strong-field-seekers can be decelerated by reversing the charging sequence so that the field is off when the molecules enter each capacitor and switched on before they leave. However, Maxwell’s equations prevent the simple focussing that weak-field seekers enjoy because it is not possible to have a static maximum of electric field in free space Wing(1)84 . This difficulty can be overcome by the method of alternating gradient focussing, whereby a sequence of alternating focussing and de-focussing lenses results in a net focussing effect. Recently a beam of strong-field-seeking CO molecules in the metastable a$`{}_{}{}^{3}\mathrm{\Pi }`$ state was slowed from 275 m/s to 260 m/s in such a decelerator Bethlem(2)02 . Here we report the deceleration and focussing of a supersonic 174YbF beam in the ground state $`X^2\mathrm{\Sigma }^+(v=0,N=0)`$ and we show that it is feasible to bring such molecules to rest. These are seven times heavier than any molecules previously decelerated, so the amount of kinetic energy to be removed is proportionally larger. If the number of decelerator stages is to remain practical each stage must remove a large amount of energy. Much more energy can be removed from heavy molecules in their ground state than in any weak-field-seeking state.

The ability to decelerate heavy molecules and ground state molecules greatly extends the range of species amenable to this powerful method of cooling and trapping. YbF is of particular interest because it offers one of the most sensitive ways to search for elementary particle physics beyond the standard model through a measurement of the electric dipole moment of the electron Hudson(1)02 . This important quantity can be measured more sensitively by slowing the molecules. Furthermore, deceleration would enhance experiments that study possible connections between the homochirality of biomolecules and parity violation Daussy(1)99 ; Ziskind(1)02  and would improve the sensitivity of nuclear anapole measurements DeMille(1)03 .

Our beam line is shown in Fig. 1(a). The supersonic source uses 4 bar of Ar, Kr or Xe mixed with a small fraction of SF6 as the carrier gas. A solenoid valve releases this gas every 0.2 s in 100 $`\mu `$s pulses that expand through a nozzle into a vacuum of $`5\times 10^5`$ mbar. A disk of Yb is positioned immediately outside the nozzle. The edge of this disk is ablated by the unfocussed 1064 nm output of a 10 ns, 20 mJ pulsed Nd:YAG laser. The Yb and SF6 react to form YbF, which becomes entrained in the gas pulse and thermalizes with it. The molecules then pass through a 1 mm-diameter skimmer, 82 mm from the nozzle, into a high vacuum. This source produces YbF beams as slow as 290 m/s with translational and rotational temperatures as low as 1.4 K Tarbutt(1)02 . A 552 nm laser beam, the velocity selector, crosses the molecular beam 108 mm downstream from the skimmer. This laser beam pumps molecules out of the ground state except during an adjustable period of typically 10-20 $`\mu `$s when it is switched off. The moment at which a molecule crosses the laser beam is mainly determined by its speed. Consequently, the resulting pulse of ground state molecules is a velocity group whose width and center are determined by the timing of the laser pulse. A further 4 mm downstream the molecules enter the decelerator. This has 12 stages, each comprising a pair of 20 mm-long stainless steel rods parallel to the beam axis z and rounded to have hemispherical ends. The rods are 6 mm in diameter, with a 2 mm gap between them, as indicated on the left of fig. 1(b). The stages are separated from each other by 10 mm gaps. They are arranged in four groups of three, alternating between vertical ($`y`$) and horizontal ($`x`$) orientation, as illustrated in fig. 1(a). Finally, a 552 nm probe laser, 258 mm from the end of the decelerator and perpendicular to the molecular beam, selectively excites one hyperfine sub-level of the ground state. The time-evolution of the fluorescence is detected by a photomultiplier.

The decelerator electrodes are operated at $`\pm 10`$ kV, giving an electric field on the axis of nearly 100 kV/cm at the center of each stage. Fig. 1(b) shows the resulting Stark shift of the molecules in a vertical stage versus transverse position $`x`$ or $`y`$. This potential focusses the molecules toward the beam axis in the $`x`$-direction and defocusses then along $`y`$. It is approximately harmonic over most of the region between the electrodes. The transverse motion can therefore be described using transfer matrices (as in optics), and the transverse oscillations are stable provided $`2<`$Tr(M)$`<2`$, where M is the transfer matrix for one unit of the periodic structure Bethlem(2)02 . We have chosen to alternate the focussing and defocussing axes after every third lens, rather than between each lens. In this way, we are able to maximize the number of deceleration stages per unit length and still satisfy the alternating gradient stability criterion.

Fig. 1(c) shows the variation of the Stark shift along the beam axis through one stage. As the beam pulse passes through the stage the voltages are switched on and off by fast high-voltage switches under the control of a 50 MHz pattern generator. Fig. 1(c) illustrates where the field might be switched on and off to decelerate a molecule moving to the right. The switching sequence that decelerates a pulse of molecules through many stages is designed so that a particular molecule, known as the synchronous molecule, reaches the same positions within each stage at the moments of switching on and off. The switch-off position is described by a phase angle $`\varphi `$ defined through the energy loss of the synchronous molecule per stage $`W_0\mathrm{sin}\varphi `$, where $`W_0`$ is the maximum possible energy loss per stage. By choosing $`\varphi <90^{}`$ we ensure that those molecules running ahead of the synchronous molecule are decelerated more, while those lagging behind are decelerated less. This phase stability is important for a decelerator because it allows a finite velocity group of molecules to propagate as a bunch Bethlem(2)00 .

The dotted curve in fig. 2(a) shows the YbF time-of-flight distribution with the decelerator off, using Xe carrier gas pre-cooled to 210 K. We deduce from this distribution that the YbF pulse has a mean velocity of 298 m/s distributed with a full width at half maximum (FWHM) of 26 m/s, corresponding to a temperature of 2.7 K. When the decelerator is turned on (solid line) at a phase angle of $`60^{}`$, approximately half the molecules gather into a narrow bunch whose speed is 287 m/s with 13 m/s FWHM. Most of the molecules inside the bunch have been decelerated, although some of the slowest molecules from the original distribution are accelerated into the bunch, as indicated by the reduction in the number of molecules at long arrival times when the decelerator is turned on. The bunch will stay together as more stages of deceleration are added to bring it to rest.

To test our understanding of the decelerator, we simulate the motion of molecules through the apparatus by integrating the axial equation of motion using the calculated potential on the axis of the beam line, which is plotted in Fig. 1(c). For the initial velocity distribution in the simulation, we take a beam whose temperature and central velocity have the values derived from the dotted data in fig. 2(a). The result of this calculation, plotted in fig. 2(b), is qualitatively very similar to the experimental result, indicating that the motion of the molecules is indeed roughly separable into axial and transverse components and showing that we understand the basic deceleration process.

Fig. 2(c) shows the results obtained with the velocity selector turned on and timed to select the central $`20\mu s`$ section of the pulse. The resulting velocity distribution remains centered on 298 m/s but now has a FWHM of only 8.5 m/s. When the decelerator is turned on, this entire distribution is shifted down in velocity by 10 m/s and the time-of-flight profile is narrowed to 70% of its original width because the decelerator slows down the faster molecules in the bunch more than the slower ones. Fig. 2(d) confirms that these results are also well described by a simple simulation of the axial motion.

We have further tested the decelerator using four other initial speeds, chosen by varying the mass and temperature of the carrier gas, with the results shown in fig. 3. In each experiment the velocity selector was used to narrow the velocity distribution to a FWHM in the range 9-14 m/s. This longitudinal velocity range lies within the phase stability range of the decelerator, meaning that the entire velocity-selected pulse can be decelerated. In all four experiments, the phase angle is set to $`60^{}`$, corresponding to a total energy loss of $`46`$ cm-1. The measured energy loss confirms this, showing that the decelerator works correctly over a wide range of initial speeds. All the beams become 20-30% narrower in their time-of-flight profile when the decelerator is switched on as a result of the phase stability. The slowest beam we have produced is shown in Fig. 3(d). Here the speed of the molecules entering the decelerator is only 286.5 m/s and the decelerator reduces their energy by 7%.

The final speed of the molecules can be adjusted by making an appropriate choice of phase $`\varphi `$ in the high voltage switching pattern, as shown by the data points in figure 4. Here we have simply taken the peak of the decelerated pulse as a measure of the speed. The solid line shows the expected change for a synchronous molecule, corresponding to an energy loss of $`W_0\mathrm{sin}\varphi `$ per stage. Crosses show the result of a numerical simulation that takes into account the entire ensemble of molecules. The uncertainty in these simulations is the same as the error bars on the data. As the phase is increased toward $`90^{}`$, the measured change of speed shows reasonable agreement with the theories. Figures 2, 3 and 4 show that this alternating gradient decelerator works well over a wide range of operating parameters and that its deceleration and phase stability properties are well understood.

Transverse confinement is important for the successful operation of the decelerator. In order to study this experimentally we measured the number of molecules reaching the detector for a range of electrode voltages and lens lengths. Although the physical length of the electrodes is fixed, we can vary the length of time for which each lens acts on the molecules by changing the turn-on point in the timing sequence (fig. 1(c)), keeping the turn-off point fixed. Measurements on a 590 m/s beam are shown in Fig. 5, where we plot the total number of molecules transmitted through the decelerator, normalized to the number when the decelerator is off for various lens lengths (5(a)) and various electrode voltages (5(b)). The dotted lines show the expected transmission of molecules through our decelerator in the absence of edge effects, lens aberrations and misalignments. This was calculated using transfer matrices, sampling over a realistic range of initial transverse velocities and positions. At low lens power, the calculated transmission dips below unity because the defocussing lenses are 10% stronger than the focussing ones. At higher power, alternating gradient focussing increases the transmission. Our data follow the expected trend, whereas the only previous study of alternating gradient deceleration of molecules showed a 20-fold shortfall in the transmission, a problem attributed to misalignments of the electrodes Bethlem(2)02 . By varying the phase angle $`\varphi `$ we can control the fraction of time that the molecules spend in the fringe fields of the lenses and we find that this alters the focussing. We believe that the fringe fields, not accounted for in the matrix model, are the cause of the small discrepancy between theory and experiment observed in fig. 5(a). This observation underlines the importance of keeping the fringe fields under review while designing alternating gradient decelerators.

We have demonstrated that heavy ground state molecules can be slowed down by using an alternating gradient decelerator. The technique we have demonstrated can be applied to a wide range of heavy polar molecules. We have used a 12-stage decelerator operated at a field of 100 kV/cm with a deceleration phase angle of $`60^{}`$ to remove 7% of the kinetic energy of YbF molecules. We were able to raise the electric field to 140 kV/cm without electric field breakdown. At this higher field, 110 stages would bring the molecules to rest. They could then be coupled into an electrodynamic trap Shimizu(1)92 ; Peik(1)99  or a wire trap Sekatskii(1)96 , or they could be switched into a weak-field seeking state and stored in an electrostatic trap e.g. a quadrupole trap Wing(1)80  or chain-link trap Shafer-Ray(1)03 . Trapping heavy molecules will provide a spectacular increase in the spectroscopic resolution available, leading to new applications in metrology, quantum chemistry and fundamental physics.

We are indebted to Henrik Haak and André van Roij for mechanical design work and to Victor Ezhov for valuable discussions. This work was supported in the UK by PPARC, EPSRC and JIF, by the EU ‘Cold Molecules’ network, and in Russia by RFBR grant 02-02-17090. HLB acknowledges a Talent fellowship from the NWO.