text:
# JETS AT ALL SCALES

## 1 Introduction

Jets are “pipelines” through which energy and matter originating from the central region of AGNs can flow out from the galaxy and, as in FRII radio galaxies, reach the radio-lobes (e.g. Begelman, Blandford & Rees 1984). The comprehension of such structures represents a great challenge for Astrophysics: the ultimate goals include the understanding of the physical processes able to produce such energetics and collimated outlows and the complex dynamics of the jet and its interaction with the environment.

In the last decade, thanks to new instruments (in particular the $`\gamma `$-ray telescope EGRET, the TeV Cherenkov telescopes and Chandra), this field has received a great impulse. In the following I briefly present our current view of the principal physical processes acting in jets, from the innermost portion of the jet close to the central “engine” to the external regions. Finally I discuss some works in progress regarding the connection between the inner portion of the jet and the outer regions.

## 2 The inner jet: Blazars

It is now generally accepted that the blazar ”phenomenon” (highly polarised and rapidly variable radio/optical continuum) is due to a relativistic jet pointing close to the line of sight (e.g. Urry & Padovani 1995). Relativistic effects strongly amplify the emitted non-thermal continuum produced close to the central BH (at $`d<0.1`$ pc, Fig.1) that can dominate over the thermal contribution from the disk or the Broad Line Region. Therefore the study of these sources offers us the unique possibility to probe the jet in its initial portion.

### 2.1 Spectral properties

The typical Spectral Energy Distribution of blazars presents two broad humps, the first peaking in the IR through X-ray band, the second in the $`\gamma `$-ray domain, extending in some sources to the TeV band (Fig.1).

The high degree of polarization clearly indicate that the low energy component is produced through synchrotron emission. The second peak is likely produced through IC scattering of soft photons by the same electron population responsible for the synchrotron emission (although other mechanisms have been proposed, see e.g. the review in Sikora & Madejski 2001). In the Synchrotron-Self Compton model it is assumed that the soft photon energy density is dominated by the synchrotron photons themselves, while the so-called External Compton models assume that soft photons coming from the external environment (disk, BLR) dominate over other possible contributions.

Although blazars form a rather inhomogeneous class from the point of view of the optical classification, the SEDs present a notable trend with the luminosity, the so-called blazars sequence (Fossati et al. 1998; Donato et al. 2001). The sources with high luminosity output ($`L10^{48}`$ erg cm-2 s-1, sources mainly classified as Flat Spectrum Radio Quasars) have both peaks located at low energies, the first in the IR region and the high energy peak around 1 MeV. On the contrary, at low power the sources (BL Lac objects, showing weak or even absent thermal features) display very high peak frequencies: in HBL, which appear to be the less powerful objects (with typical luminosities $`10^{42}`$ erg cm-2 s-1) the first peak falls in the UV-soft X-ray band and the high energy component can peak near the TeV region.

This behaviour has been interpreted (Ghisellini et al. 1998, recently revisited in Ghisellini et al. 2002) as the result of a balance between a acceleration an cooling suffered by the electrons producing the radiation at the peaks. In the assumption that the acceleration rate is constant in all the sources, the sequence can be successfully interpreted as due to an increasing cooling rate (dominated in most of the sources by IC losses) from the low-power BL Lacs to the powerful FSRQs.

### 2.2 Jet power and matter content

SEDs can be satisfactorily reproduced by simple emission models (e.g. Ghisellini et al. 1998, Tavecchio et al. 1998), yielding the value of the main physical quantities of the jet, such as electron density and energy, magnetic field, size of the region. In particular measuring the X-ray spectra and adapting a broad band model to their SEDs yields reliable estimates of the total number of relativistic particles involved, which is dominated by those at the lowest energies. This is interesting in view of a determination of the total energy flux along the jet (e.g. Celotti et al. 1997, Sikora et al. 1997). The total ” kinetic” power of the jet can be written as:

$$P_{\mathrm{jet}}=\pi R^2\beta cU\mathrm{\Gamma }^2$$
(1)
where $`R`$ is the jet radius, $`\mathrm{\Gamma }`$ is the bulk Lorentz factor and $`U`$ is the total energy density in the jet, including radiation, magnetic field, relativistic particles and eventually protons. If one assumes that there is 1 (cold) proton per relativistic electron, the proton contribution is usually dominant.

In Fig. 2 the derived radiative luminosity $`L_{\mathrm{jet}}`$ and kinetic power of the jet $`P_{\mathrm{jet}}`$ for a group of sources with sufficiently good spectral information are compared. The ratio between these two quantities gives directly the “radiative efficiency” of the jet, which turns out to be $`\eta 0.1`$, though with large scatter. The line traces the result of a least-squares fit: we found a slope $`1.3`$, suggesting a decrease of the radiative efficiency with decreasing power.

### 2.3 The accretion-jet connection

Two main classes of models for the production of jets has been proposed. The first class considers the extraction of rotational energy from the black hole itself (Blandford & Znajek 1977), while the in the models of the other class one considers MHD winds associated with the inner regions of accretion disks (Blandford & Payne 1982). In both cases the extractable power critically depends on the value of the magnetic field close to the BH horizon, which is determined by the rate at which matter is falling onto the BH. Therefore a connection between the disk luminosity (proportional to $`\dot{M}`$) and the jet power (proportional to $`B\dot{M}`$) is naturally expected.

A simple example can be discussed for the case of BZ model and spherical accretion. The power extractable in the BZ power can be written as:

$$P_{BZ}\frac{1}{128}B_0^2r_g^2a^2c$$
(2)
where $`r_g`$ is of the order of the gravitational radius of the BH. Assuming maximal rotation for the black hole ($`a=1`$), the critical problem is the estimate of the intensity reached by the magnetic field threading the event horizon, which must be provided by the accreting matter. Using a spherical free fall approximation with $`B_0^2/8\pi \rho c^2`$ one can write:

$$P_{BZ}g\dot{M}c^2$$
(3)
where $`P_{acc}=\dot{M}c^2`$ is the accretion power and $`g`$ is of order 1/64 in the spherical case, but in fact it is a highly uncertain number since it also depends on the field configuration.

Estimates of the power of jets and of the associated accretion flows can therefore be crucial to shed light on the jet-disk connection . In a pioneering work Rawling & Saunders (1991) addressed this question studying a large sample of radio galaxies, finding a good correlation between the two. Celotti, Padovani & Ghisellini (1997) first investigated the relation between jet and disk using the direct VLBI radio emission of the jet close to the nucleus for the jet power estimation. They found a suggestive hint of correlation between these two quantities, although the statistical significance was too low to draw a firm conclusion.

For most of the sources considered in the previous section we have also information on the amount of radiation released by the accretion processes, estimated either through the observed emission lines or, in few cases, from the observed blue bumps, direct signature of the disk. For BL Lac sources we can only infer upper limits on the disk luminosity, but for the case of BL Lac itself, which shows a bright line. In Fig. 2 we compare the luminosity of the jet, $`L_{\mathrm{jet}}`$, which is a lower limit to $`P_{\mathrm{jet}}`$, with the luminosity of the disk, $`L_{\mathrm{disk}}`$. A first important result is that on average the minimal power transported by the jet is of the same order as the luminosity released in the accretion disk. More precisely the data we have used suggest $`L_{jet}L_{disk}`$ at high luminosities and $`L_{jet}>L_{disk}`$ at intermediate and low luminosities.

Assuming that $`\eta =10^1`$ as estimated above and that the accretion takes place with the standard radiative efficiency $`ϵ0.1`$, the near equality of $`L_{jet}`$ and $`L_{disk}`$ then requires $`P_{\mathrm{jet}}=P_{\mathrm{acc}}`$.

On the other hand, a dominance of $`L_{jet}`$ over $`L_{disk}`$ at lower luminosities could be attributed to a lower value of $`ϵ<<0.1`$ which may be expected if the accretion rate is largely sub-Eddington. In the latter case the range in luminosities spanned by Fig. 2 should be mainly a range in accretion rates rather than a range in black hole masses. For instance the minimum jet powers of three of the BL Lacs in Fig 2 are around $`10^{44}`$ erg/s which suggests $`P_{jet}10^{45}`$ requiring a mass of $`10^7`$ for critical accretion rate. Since the disk luminosity is less than $`10^{42}`$, if the accretion rate is 1% Eddington the implied mass is again $`10^9`$ $`M_{}`$.

## 3 Extended X-ray jets: Chandra results

Although very common in the radio band, before the launch of the Chandra satellite in 1999 only a handful of extragalactic kpc scale jets were known to emit X-rays. Among them the bright and prominent jets in 3C273 , M87, Cen A, studied with EINSTEIN and ROSAT. With the superior sensitivity and, especially, spatial resolution of Chandra numerous jets have been detected, triggering a new intense theoretical and observational work. Even in the source selected for the first light, the distant ($`z=0.6`$) quasar PKS 0637-752, a prominent jet has been discovered (Chartas et al. 2000).

Soon after the discovery of the X-ray jet in PKS 0637-0752 other jets have been detected both in radio galaxies and in quasars. A recent census (see e.g. the WEB site maintened by D. Harrisaaahttp://hea-www.harvard.edu/XJET/) reported 25 jets detected in X-rays. Different classes of AGNs are represented in the sample: the most numerous group is that of radio-galaxies (both FRI and FRII), but a large fraction is composed by powerful radio-loud QSOs. Most of the jets have also an optical counterpart, usually detected by HST.

The first problem posed by observations is the identification of the emission mechanism responsible for the production of X-rays. Since these jets are known to emit in radio, the first candidate mechanism is synchrotron emission. In some cases (in particular for jets in radio-galaxies) this interpretation is consistent with the data, but in other several cases (especially in quasars) it fails to explain the observational evidence. This is the case of the first jet discovered in PKS 0637-0752. For this specific case Tavecchio et al. (2000) and Celotti, Ghisellini & Chiaberge (2000) proposed that the emission responsible for the observed radiation is the IC scattering of the CMB radiation by relativistic electrons in the jet.

In order to study in a more systematic way the properties of the X-ray jets we proposed a series of combined Chandra-HST observations of a selected group of quasar known to have radio jets. We also started a program to image these jets uniformly at radio wavelengths using the VLA, VLBA, and MERLIN. The work is still in progress. Details for a first group of selected sources can be found in Sambruna et al. (2001). An example is reported in Fig. 3.

## 4 From the Blazar region to the kpc scale

The possibility to constrain the physical state of the plasma in the jets both at Blazars and kpc scale could offer us the interesting opportunity to shed some light on the evolution of the jet from very small scales, close to the central engine, to the outer regions, where the jet is starting to significantly decelerate. In fact we are starting to do this work for a small subset of the Chandra jets for which good data for both regions are available (Tavecchio et al., in prep).

A first important result is that the powers independently derived for both regions are very close, supporting the overall metodology and suggesting the idea that the jet does not substantially dissipate its power until its end. Another interesting suggestion derives considering the internal pressure of the knots derived with our modelling. The derived values are in fact close to $`10^{11}`$, suggestively close to the pressure of typical hot gas halos found around FRIs (Worrall et al. 2000).

The results discussed above seems to point toward a simple scenario for the dynamical evolution of the jet from inner to the outer regions. Clearly the pressure in the blazar jet is larger than any possible confining gas. This strongly suggests that, after the dissipation region, the jet expands freely until the internal pressure (gas and/or magnetic) reaches the pressure of the (supposed) hot gas of the galactic halo. At this point a reconfinement shock forms (the jet is highly supersonic), particle are accelerated, and the jet starts to decelerate (Sanders 1983).

## References


format: text

text:
# A New Mechanism of Model Membrane Fusion Determined from Monte Carlo Simulation

## I  Introduction

Although membrane fusion is a fundamental biological process of importance in fertilization, synaptic release, intracellular traffic, and viral infection, its basic mechanism is not well understood. Much of the literature has focused on fusion proteins whose function is, inter alia, to overcome the energetic cost of bringing the bilayers to be fused to within a small distance of one another, a step which places the membranes under tension Chen and Scheller, (2001). There is accumulating evidence, however, that the subsequent stages in the fusion pathway, the interruption of the integrity of the bilayers and the molecular rearrangements that lead to the formation of the fusion pore itself, are essentially lipidic in nature (Zimmerberg and Chernomordik,, 1999; Lentz et al.,, 2000). A consequence of this view is that the fusion process can be studied, both experimentally and theoretically, utilizing simple model membrane systems. Knowledge of the fusion mechanism in these simpler systems would illuminate additional roles that the proteins need to play in biological fusion.

The theoretical treatment of membrane fusion has, almost without exception, been restricted to phenomenological models which decribe the bilayer, not in terms of the microscopic architecture of its components, but rather in terms of the macroscopic elastic properties of its monolayers. The common assumption is that these elastic moduli are uniform and independent of membrane deformations Safran, (1994). Although attractive mathematically, this approach has its limitations. For instance, it is not clear whether the expansion of the membrane free energy to second order in deformations is sufficient to describe the highly curved intermediate structures, which may be involved in fusion. Additional approximations must be introduced to calculate the properties of junctions of bilayers, which are not well described by simple bending deformations. The energy of these structures has proven to be particularly sensitive to the approximation used in their description Siegel, (1993); Kuzmin et al., (2001); Kozlovsky and Kozlov, (2002); Markin and Albanesi, (2002). Importantly, application of these approaches requires one to assume a particular fusion pathway. The only pathways considered to date has been limited to variants of one hypothesis Markin and Kozlov, (1983); Kozlov and Markin, (1983); Chernomordik et al., (1985). One starts with two bilayers in close apposition. Lipids in the facing, or cis, layers rearrange locally and bridge the aqueous gap between the bilayers. This results in the formation of an axially symmetric stalk. In most versions, the stalk then expands radially and the cis layers recede. The trans layers make contact and produce an axially symmetric hemifusion diaphragm. Nucleation of a hole in this diaphragm completes the formation of an axially symmetric fusion pore. Because of the evolution of the stalk into a hemifusion diaphragm in this model, we shall refer to it as the “hemifusion mechanism.” Because only variants of the hemifusion mechanism have been examined, and because the theory is phenomenological, one does not know a priori in what systems this pathway may be the most favored, or under what conditions. Some insight is gained by comparison with experiment which shows this hypothesized mechanism to be consistent with a wide range of experimental observations of biological lipids Monck and Fernandez, (1996); Zimmerberg and Chernomordik, (1999); Jahn and Grubmüller, (2002). However there is no direct evidence to confirm that this particular pathway is that taken either by biological or laboratory-prepared model membranes.

In light of the above, it would certainly be desirable to examine the fusion pathway in a system whose components are described by a microscopic model. Such examination has begun recently. A minimal model, consisting of rigid amphiphiles of one hydrophilic and two hydrophobic segments, and with no explicit solvent, was studied with Brownian dynamics simulations Noguchi and Takasu, (2001). At the same time, a model of more complex, flexible, chain molecules, widely employed in the polymer community, was used by us to study bilayers composed of amphiphilic, diblock copolymers in a hydrophilic solvent Müller et al., (2002). Such copolymers are known, in fact, to form bilayer vesicles which can undergo fusion Discher et al., (1999). This system was studied by means of Monte Carlo simulations. Both theoretical studies observed the formation of the initial stalk, but found that the subsequent fusion pathway was not the usual hemifusion mechanism, but involved intermediates that broke the axial symmetry. In particular, off to the side of the initial stalk, the formation of small pores in each of the fusing bilayers was clearly seen. (We shall refer to these small pores which span one bilayer only as “holes” to avoid confusion with fusion pores which span both bilayers.) It is intriguing that the two studies observed the same fusion pathway even though the architecture of the constituents of the two systems differed considerably, sharing little other than the generic property of being amphiphilic and capable of bilayer self assembly.

The two investigations gave a first glimpse of a fusion pathway which differs from the hemifusion mechanism, but did not provide a great deal of quantitative detail. In this paper we present an extensive study of the same microscopic model we employed previously, and offer sufficient quantitative evidence to substantiate our earlier observations.

Naturally we are concerned with the question of whether the fusion pathway we observe in our model system is relevant to membrane fusion in biological systems. The architecture of the components in our system obviously differs greatly from those of biological lipids, and it is not clear how one should compare the systems. We make such an attempt by calculating several dimensionless ratios which can be formed from membrane parameters and comparing those in our system with ratios characteristic of vesicles formed of block copolymers, and of liposomes. (See Table I below.) Ultimately we cannot be sure of the systems to which our results apply and under what conditions, save the very particular ones that we have simulated for the particular case of block copolymers. In this sense, our results must be evaluated in the same way as those from the phenomenological theories; they must be compared to experiment. We do so in the Discussion. In particular we note that our mechanism predicts that the fusion rate depends on lipid architecture and membrane tension, that there is mixing of lipids in the cis leaves before mixing of contents, and that there is also mixing of lipids between cis and trans layers. Of most interest, our mechanism predicts that transient leakage is causally linked to the process of membrane fusion.

## II  Simulation details

Simulation of membrane fusion in a fully chemically realistic model would be most valuable, because it could provide information about specific structural changes on the atomic level. This would be particularly important if changes in molecular conformations entailed a qualitative spatial redistribution of hydrophilic and hydrophobic segments. Unfortunately, the simulation of atomistically faithful models can only follow the time evolution of a few hundreds of lipid molecules over a few nanoseconds even on state-of-the-art supercomputers. Given that the time scale of membrane fusion is on the order of milliseconds and involves lengths on the order of a few tens of nanometers, an atomistic simulation of the fusion process is not yet feasible and one has to resort to coarse-grained models.

Coarse-grained models of amphiphilic chain molecules have been used with great success to investigate common features of self-assembly. Such models retain only those molecular properties that are necessary for self-assembly, such as the connectivity of hydrophilic and hydrophobic portions along the amphiphilic molecule, and the mutual repulsion between these different kinds of segments, and ignore specific chemical or electrostatic interactions. The usefulness of this approach rests on the observation that chemically very different systems, such as biological lipids in aqueous environment and block copolymers in a homopolymer environment, exhibit a common phase behavior and similar structural patterns on length scales comparable or larger, than the molecular size. The self-assembly of amphiphiles into bilayer membranes itself is an example of a universal behavior, i.e. one which does not depend on fine details of the underlying architecture. It has been successfully studied by coarse-grained models Shillcock and Lipowsky, (2002). We expect that all membranes can be caused to fuse, however there may be several different pathways which are taken by different systems under different conditions. Our purpose here is to demonstrate one path which is taken in a system modeled microscopically.

We employ the bond fluctuation model Carmesin and Kremer, (1988) of a polymer chain, which has been used previously to study pore nucleation in a symmetric bilayer membrane under tension Müller and Schick, (1996). Much is known about the structure and thermodynamics of this model, and the parameters can be mapped onto the standard Gaussian chain model of a dense mixture of extended molecules. In this three-dimensional lattice model, each segment occupies a lattice cube. No two occupied cubes can share any corner, a rule that mimics hard-core repulsion interaction. Furthermore this ensures that the lattice spacing is sufficiently smaller than the width of interfaces so that the effect of the discretization of space is minimal. To ensure that the chain of segments can not intersect itself, the segments are connected by bond lengths that cannot be too large. In particular, neighboring segments along the chain can be connected by one of 108 bond vectors of lengths $`2,\sqrt{5},\sqrt{6},3`$ or $`\sqrt{10}`$ measured in units of the lattice spacing $`u`$. The angles between adjacent monomers can take on any of 87 values. The large number of bond vectors and the extended segment shape allow a rather faithful approximation of continuous space, while retaining the computational advantages of lattice models. The amphiphilic molecules consist of $`11`$ hydrophilic segments and $`21`$ hydrophobic segments. This asymmetry mimics the ratio of head and tail size in biologically relevant lipid molecules, and is slightly smaller than employed by us previously Müller et al., (2002). We reduced, in this study, the asymmetry of the molecules so that a solvent-free system not only would be in a lamellar phase ($`L_\alpha `$), but would also be further than in our previous study from the boundary separating the lamellar and inverted-hexagonal ($`H_{II}`$) phases. The solvent in our system is represented by chains of $`32`$ hydrophilic segments, i.e. we conceive a hydrophilic chain as a small cluster of solvent molecules, just as in other coarse-grained modeling Shelley et al., (2001) The mean head-to-tail distance of the amphiphiles and solvent molecules is $`17u`$. Like segments attract each other and unlike segments repel each other via a square well potential which comprises the nearest $`54`$ lattice sites. Each contact changes the energy by an amount $`ϵ=0.17689k_BT`$. This corresponds to an intermediate segregation $`\chi N30`$ in terms of the Flory–Huggins parameter $`\chi `$. If we increased the incompatibility much more, we would reduce the interfacial width between hydrophilic and hydrophobic segments to the order of the lattice spacing and the local structure of the lattice model would become important. If we decreased the incompatibility, we would reduce the clear segregation between hydrophilic and hydrophobic regions. Similarly if we replaced the solvent homopolymers by monomers, we would effectively reduce the incompatibility Matsen, (1995), and again reduce the segregation between the diblock and solvent hydrophilic segments and the diblock hydrophobic segments. Were we to increase the incompatibility to restore the desired degree of segregation, we would again reduce the interfacial width of the membrane to an extent that lattice effects would become important.

Monte Carlo simulations are performed in the canonical ensemble, except for some runs described in section III. The segment number density, i.e. the fraction of lattice cubes occupied by segments, is fixed at $`\rho =1/16`$. The conformations are updated by local segments displacements and slithering-snake-like movements. The different moves are applied with a ratio $`1:3`$. We count one attempted local displacement per segment and three slithering-snake-like attempts per molecule as four Monte Carlo steps (MCS). This scheme relaxes the molecular conformation rather efficiently. The latter moves do not mimic a realistic dynamics of lipid molecules and we cannot identify straightforwardly the number of Monte Carlo steps with time. The density of hydrophilic and hydrophobic segments, however, is conserved so that the molecules move diffusively. Moreover, the molecules cannot cross each other during their diffusive motion. In this sense we have a slightly more realistic time evolution on local length scales than in dissipative particle dynamics simulations Shillcock and Lipowsky, (2002), but Monte Carlo simulations cannot include hydrodynamic flows, which might become important on large length scales. At any rate, we do not expect the time sequence to differ qualitatively from that of a simulation with a more realistic dynamics on time scales much larger than a single Monte Carlo step. Most importantly, fusion is thought to be an activated process, therefore the details of the dynamics only set the absolute time scale, but the rate of fusion is dominated by free energy barriers encountered along the fusion pathway, which are independent of the actual dynamics used.

## III  Preparation and properties of a single bilayer

It seems clear that bilayers that are under no stress will not undergo fusion, as there is no free energy to be gained by doing so. So to promote fusion, we have subjected the studied bilayers to lateral tension. This has been done by providing the system with fewer molecules than are needed to span the given area of our sample cell with bilayers that are tensionless. Of course we need to know just how many molecules are needed to make a tensionless bilayer that spans the cell. To determine properties of the tensionless bilayer, we made use of the definition of the tension in this liquid-like bilayer as the derivative of the free energy with respect to the bilayer area at constant temperature and particle number. We therefore investigated an isolated bilayer with a straight, free edge. A simulation cell of size $`64u\times 200u\times 64u`$ with periodic boundary conditions in all dimensions was used. The bilayer, oriented in the $`xy`$ plane, spanned the system in the short, $`x`$, direction, but did not span the system in the long, $`y`$, direction. Its extension in this direction adjusted itself until it neither grew nor shrank. Thus the surface tension, $`\gamma `$, of the bilayer was zero. This vanishing value includes, of course, the contributions from the fluctuations of the bilayer. Even though the tension vanished, these fluctuations of the mid-plane were not very large due to the stiffness of the rather small patch of membrane considered. A typical snapshot of the bilayer configuration is shown in Fig. 1(a). Rearrangement of amphiphiles at the bilayer free edges is clear. The average profile along the $`y`$ axis, the long axis of the bilayer, is presented in panel (b). To obtain it, we have averaged the profiles along the $`x`$ and $`z`$ direction and estimated the instantaneous angle the bilayer makes with the $`z`$ direction (to correct for the difference between projected and true area). We observe for these laterally averaged profiles that the edge of the bilayer is slightly thicker than the middle, the increase is about $`7\%`$ for the amphiphilic segment density and about $`16\%`$ for the density of hydrophobic segments. Away from the edge, the densities decay exponentially to those of the uniform bilayer (i.e without an edge), and we estimate the thickness of the tensionless bilayer from that in the middle, finding it to be $`d_0=31u`$.

The profiles across a single bilayer of thickness $`d_0=31`$ are shown in Fig. 1(c). They were obtained by simulation in a cell $`40u\times 40u\times 80u`$ in which the bilayer spanned both short directions. One sees that hydrophobic and hydrophilic regions are clearly separated, but there is some interdigitation of the hydrophobic tails emerging from the opposing monolayers.

Knowing the thickness of the tensionless bilayer, we know the number of molecules needed to span the simulation cell with such a bilayer, and can control tension by varying the number of molecules introduced into the cell. We cannot determine this tension, as one might in a Molecular Dynamics simulation, from the excess tangential stress in the interfacial zone because we employ a lattice model. Nevertheless, we can determine the tension purely from thermodynamic relations. To do so, we assembled a single bilayer in a system of size $`156u\times 156u\times 64u`$, where the bilayer spanned the system in the $`xy`$ plane. Using semi-grand canonical identity switches between amphiphiles and solvent, we controlled the exchange potential $`\mathrm{\Delta }\mu `$ between the species and monitored the thickness of the bilayer (measured by the areal density of amphiphiles). The dependence of the bilayer tension on the chemical potentials of the amphiphile, $`\mu _C`$, and solvent, $`\mu _A`$, is given by the Gibbs absorption isotherm Davis, (1996):

$$L^2\mathrm{d}\gamma =\delta n_C\mathrm{d}\mu _C\delta n_A\mathrm{d}\mu _A\delta n_C\mathrm{d}\mathrm{\Delta }\mu $$
(1)
where $`\delta n_C`$ and $`\delta n_A`$ are the excess number of molecules in the bilayer. In the last approximation we have assumed that the liquid is incompressible i.e. $`\delta n_A\delta n_C`$, and the solubility of the amphiphile in the solvent is vanishingly small. Results of the simulation for the number of amphiphiles $`\delta n_C`$ as a function of the exchange potential $`\mathrm{\Delta }\mu =\mu _C\mu _A`$ are shown in Fig. 1(d). Using the thickness of the tensionless bilayer, we can estimate the tension of an arbitrary bilayer as a function of exchange potential or of thickness by integrating Eq. (1). The results, in reduced units of bare A-B homopolymer interfacial tension $`\gamma _0=0.068k_BT/u^2`$, are shown in the inset of Fig. 1(d). Dashed lines in Fig. 1(d) and the inset correspond to the tensionless bilayer. Comparison of the relevant structural and elastic properties of the polymersomes, liposomes and simulated membranes is provided in Table 1.

We are now in a position to simulate bilayers under a given tension in the canonical ensemble. Knowing the area of our simulation cell, and the segment density, we add the number of amphiphiles which will produce a bilayer of a given thickness. From Fig. 1(d), we know what tension is placed on this bilayer. For our study of two bilayers under tension, we have chosen their thickness to be $`d=25u`$, smaller than the thickness $`d_0=31u`$ of the tensionless bilayer. This corresponds to a tension of the order of $`\gamma /\gamma _00.75`$ and an area expansion, $`\mathrm{\Delta }A/A_00.19.`$ We know from our simulations that a single bilayer of the thickness chosen, $`d=25u`$, is metastable on the time scale of fusion, i.e. the small holes, which appear transiently, do not grow past their critical size on the time scale of fusion in our simulations.

## IV  Preparation of a System of Two Bilayers

We begin with a system containing only amphiphiles. It is $`156u\times 156u\times 25u`$ with periodic boundary conditions in the long directions, and hard, impenetrable, walls in the short direction. They attract the hydrophilic portion of the amphiphile and repel the hydrophobic portion. These interactions extend over two layers nearest to the wall, and each contact changes the energy by $`0.6k_BT`$. The amphiphiles assemble into a bilayer structure which is free of defects.

Two of these flat bilayers are then stacked on top of each other with a distance of $`\mathrm{\Delta }`$ between them, and are embedded into a simulation cell with geometry $`156u\times 156u\times 126u`$. There are no walls at this point, and periodic boundary conditions are utilized in all three dimensions. The conditions of flat bilayers mimic the approach of two vesicles whose radii of curvature are much larger than the patch of membrane needed for fusion. The solvent of homopolymers is then inserted into the simulation cell via grand canonical, configurational bias, Monte Carlo moves at infinitely large chemical potential of the homopolymer until the segment number density of $`\rho =1/16`$ is reached. The initial distance $`\mathrm{\Delta }`$ between the bilayers translates into the thickness of the residual solvent layer between the two membranes. We have carried out the most extensive series of runs with $`\mathrm{\Delta }=10u`$ and unless specified otherwise, all our results are for that separation. Because the solvent homopolymers are flexible coils, and each represents a cluster of solvent molecules, many layers of solvent segments are represented between the bilayers at this separation. In our previous simulations Müller et al., (2002), we set $`\mathrm{\Delta }=0`$ and observed qualitatively similar behavior as we do with this larger separation. We increased the separation for this extensive study because, as expected, the rate of fusion events decreased (see Sec. V. below) and this allowed us to observe the sequence of structural rearrangements more clearly than in our previous work. The separation chosen, a bit less than half the thickness of one bilayer, is comparable to the separation at which fusion occurs when mediated by hemagglutinin Flint et al., (2000). A snapshot of the two bilayers is shown in Fig. 2. Hydrophobic and hydrophilic segments of amphiphiles are shown as dark and light gray spheres. For clarity solvent segments, which are present in the simulation, are not shown

Thirty-two independent starting configurations were prepared, each containing $`194,688`$ segments corresponding to about $`3,613`$ amphiphiles and $`3,708`$ solvent molecules. After every $`25,000`$ Monte Carlo steps, a configuration was stored for further analysis. Ten thousand hours of CPU time were utilized in the course of this investigation, with thirty two processors running for about two weeks.

## V  The Process of Fusion

It is straightforward to monitor the internal energy of the system during the simulation because this energy arises solely from contacts between segments, and the locations of all segments are known. (In contrast, the free energy cannot be obtained directly.) We show in Fig. 3 the behavior of the internal energy of two systems, one separated by a distance $`\mathrm{\Delta }=4u`$ (squares), and the other with $`\mathrm{\Delta }=10u`$ (circles). The energy is plotted in units of $`k_BT`$, as a function of time, in units of $`25,000`$ Monte Carlo steps. The energy initially decays, which reflects the equilibration of the system. During this initial relaxation of the starting configuration the interface between the bilayer and the solvent adjusts locally. The time scale of this initial relaxation (less than $`25,000`$MCS) is independent of the distance between the bilayers, and is about two orders of magnitude smaller than the time scale on which the fusion pore forms. Due to this separation of time scales between initial relaxation and fusion we do not expect the preparation of the starting configuration to affect the fusion process. Similarly we do not expect our results to depend on our particular choice of relaxation moves, as other choices would also lead to relaxation of the bilayers which takes place on a much shorter time scale than does fusion.

After the initial relaxation, two subsequent time regimes can be identified. First the energy rises slowly. Two mechanisms contribute to this increase of the energy. On the one hand, capillary waves of the hydrophilic/hydrophobic interfaces become thermally excited. They increase the effective interface area and thereby lead to a slow increase of the energy. Additionally, undulations result in the formation of stalks and holes. We shall discuss the details of this process in the next subsection. Later, around $`320\times 25000`$ MCS the energy decreases rapidly. This final decrease of the energy results from the fusion of the membranes which releases some of the tension stored in them. As noted above, the fusion occurs more rapidly the closer the bilayers, as expected. The increase in energy preceding fusion reflects the formation of fusion intermediates, the focus of our study.

The inset of Fig. 3 shows the fluctuations in the energy, i.e., the fluctuations between the 32 different runs at equal time. Strong fluctuations indicate energy differences between the independent runs. The peak at around $`t320`$ indicates that some systems have already formed a fusion pore (and therefore have a lower energy) while others systems have only stalks and holes (and therefore have a higher energy). The vertical bar indicates the time we have chosen to indicate on several figures the onset of fusion. The width of the peak provides an estimate for the spread of the time at which a fusion pore appears.

### V.1 The stalk and associated hole formation

During the initial stage of simulations the fluctuating bilayers collide with one another frequently and sometimes form small local interconnections. For the most part, these contacts are fleeting. Occasionally we observe sufficient rearrangement of the amphiphiles in each bilayer to form a configuration, the stalk, which connects the two bilayers, (see Fig. 4(a)) and which is not as transient. Such a stalk was hypothesized long ago to be involved in the initial stages of fusion Markin and Kozlov, (1983); Kozlov and Markin, (1983). In contrast to stable arrays of stalks which have been observed in block copolymer melts Disko et al., (1993) and in lipid systems Yang and Huang, (2002), those we see are isolated, and increase the free energy of the system. We infer the latter from two observations: that the appearance of stalks is correlated with the increase in the internal energy of the system as a function of time shown in Fig. 3; that some stalks vanish without proceeding further to a fusion pore. Thus it appears that the stalk represents a local minimum along the fusion pathway. Density profiles of the hydrophilic and hydrophobic parts of the amphiphiles in the presence of the stalk, and obtained by averaging over configurations, are shown in Fig. 5. The dimples in the membranes at each end of the stalk axis are notable. What can barely be seen is a slight thinning of each bilayer a short distance from the axis of the stalk.

After stalks are formed, the rate of formation of holes in either of the two bilayers goes up markedly. This can be seen in Fig. 6. in which we plot the fractional area of holes as a function of time for the system of two apposed bilayers, and compare it to the rate of hole formation in an isolated bilayer. In contrast to the large increase in the area of holes formed in the apposed bilayers at “time” $`t=200\times 25000`$ MCS, the fractional area in single bilayers fluctuates somewhat about an average value which is rather constant over time at a value of approximately $`0.004`$. Comparison with Fig. 3 shows that the increase in the rate of hole production in the apposed bilayers in this system with bilayer spacing $`\mathrm{\Delta }=10u`$ is correlated in time with the decrease in the energy of the system, and it is reasonable to infer that the decrease in energy is caused by the production of holes and, later, the fusion pore. Similarly, during the time before this increase in hole production, stalks are forming, and it is also reasonable to infer that the increase in energy is due to their formation.

The locations of stalks and holes are correlated; holes form close to the stalks, and the stalk elongates and moves to surround the hole. A snapshot of this is shown in Fig. 7(a) and (d). In both snapshots an elongated stalk is seen and a small hole is formed in the upper bilayer next to the stalk. The extent to which holes are, on average, found close to a stalk can be determined from the hole-stalk correlation function

$$g(r)\frac{_{r_s,r_h}\delta (|𝐫_s𝐫_h|r)P_{sh}(𝐫_s,𝐫_h)}{_{r_s,r_h}\delta (|𝐫_s𝐫_h|r)}$$
(2)
where $`P_{sh}(𝐫_s,𝐫_h)`$ is the joint probability that the lateral position $`𝐫_s`$ is part of a stalk and $`𝐫_h`$ is part of a hole, and $`\delta (r)`$ is the Dirac delta function. The value of $`g(r)`$ at large distances $`r`$ is proportional to the product of the areal fraction of holes and stalks. This correlation function is shown in Fig. 8. The scale of $`g(r)`$ increases with time indicating the simultaneous formation of stalks and holes. The figure shows that the correlation peaks at a distance of about $`16u`$, and falls rapidly at larger distances. (Recall that each bilayer has an average thickness of $`25u`$.)

It is not difficult to understand why the presence of a stalk promotes hole formation. First, if the hole forms close to a stalk, then the line tension, or energy per unit length $`\lambda `$, of that part of the hole near the stalk is significantly reduced. This can be seen from the schematic in Fig. 9. In the upper part of the figure, we show a hole which has formed far from a stalk, while in the lower, we show a hole which has formed close to one. It seems clear that the line tension in the latter is reduced simply due to the reduction of curvature of the hydrophobic-hydrophilic interface. The second reason that the stalk formation encourages the appearance of holes is due to the slight thinning of the membrane in the vicinity of the stalk to which we alluded earlier. Further it has been suggested recently that the local surface tension in the neighborhood of a defect, such as a stalk, is increased significantly Kozlovsky et al., (2002) making such a location the likely site of hole formation.

Now that one hole has formed next to the stalk, and the stalk has begun to surround it, two other events occur to complete the formation of the fusion pore. They are: i) a second hole forms in the other bilayer, ii) the stalk surrounds the hole(s) to form the rim of the fusion pore. We have observed these steps to occur in either order, and will briefly discuss them separately.

### V.2 Pathway 1. Rim formation followed by appearance of a second hole

In this scenario, a hole appears in one bilayer and the stalk completely surrounds it rather rapidly. A snapshot of the system in this configuration is shown in Fig. 7(b). This looks very much like a hemifusion diaphragm which has been suggested by many authors as an intermediate stage in fusion Markin and Kozlov, (1983); Chernomordik et al., (1985); Siegel, (1993). However, this diaphragm is quite different from the usual hemifusion one that consists of two trans monolayers of the fusing membranes. In contrast, the diaphragm we observe is made of one of the pre-existing bilayers; that is, it is made of cis and trans leaves. The appearance of a hole in this diaphragm, as shown in Fig. 7(c), and its expansion completes the formation of the fusion pore.

### V.3 Pathway 2. Appearance of second hole followed by rim formation

In this scenario, a hole appears in one bilayer and, before the stalk completely surrounds it, a second hole appears in the other bilayer. The stalk tries to surround them both, and aligns them in doing so. In Fig. 7(e) we show one stage in this process. One sees a large hole in the upper bilayer. A small hole is formed in the lower bilayer next to the stalk. Eventually, the stalk aligns and completely encircles the holes (see Fig. 7(f)) to form the final fusion pore shown in Fig. 4(b). Again, the driving force for the stalk to surround the two holes is the reduction in their (bare) line tension. Because the stalk aligns and surrounds two holes, we observe this pathway to be somewhat slower than that of pathway 1 in which the stalk need only surround one hole.

Once the fusion pore has formed, by either of the above mechanisms, it expands, driven by the reduction in surface tension. The growth of the fusion pore eventually slows and ends as the pore reaches its optimum size determined by the finite size of our cell.

## VI Discussion

We have carried out extensive Monte Carlo simulations on the fusion of two bilayer membranes comprised of amphiphilic molecules immersed in solvent. The amphiphiles and solvent are modeled by copolymers and homopolymers respectively. The membranes are under tension. The mechanism of fusion that we see begins with a stalk, as posited years ago, and incorporated in almost all fusion scenarios. However what follows after stalk formation is different from all other mechanisms which have been proposed save that presented independently by Noguchi and Takasu Noguchi and Takasu, (2001). In particular, the fusion intermediates we see break the axial symmetry which has been assumed in almost all previous calculations. We observed that the stalk destabilizes the bilayers by catalyzing the creation of small holes in them. We argued that the mechanism behind this is quite simple: the energy per unit length of the edge of a hole is reduced when the edge is adjacent to a stalk. For the same reason, the stalk will try to surround the hole formed in one bilayer once the hole has appeared. Two slightly different pathways to the final fusion pore were observed differing only on whether the hole in the second bilayer, which is necessary for complete fusion, appears before or after the stalk completely surrounds the first hole.

The question now arises as to whether the pathway we see in the model system is that which occurs in biological fusion. There are many differences between the model studied and a biological system. Perhaps the most obvious is that we have modeled flexible, single chain block copolymers, not lipids with two semiflexible tails and a rigid head. How is one to determine whether these architectural differences are significant? It is useful to recall that phenomenological theories completely ignore the architecture of the membrane constituents and encapsulate their effects in a small number of parameters which enter the theory, such as the monolayer spontaneous curvature and bending modulus. In that same spirit, we can extract from our simulation those same parameters and compare dimensionless ratios of them to those of other systems. We have done that, and presented the results in Table I. One sees that the values we obtain are reasonable. The ratio of the bilayer compressibility modulus to the hydrophilic/hydrophobic interface tension, $`\kappa _a/\gamma _0`$, closer to that of liposomes than of polymersomes. The reverse is true for the ratio of the monolayer bending modulus to the product of surface tension and the square of the hydrophobic thickness $`\kappa _b/\gamma _0d_c^2`$. One line in the table deserves comment, that for the experimental values of the bilayer area expansion, $`\mathrm{\Delta }A/A_0`$ quoted at rupture, (the “critical values”). That for liposomes is smaller than that for polymersomes at rupture, which is equal to the bilayer area expansion we utilized. However the values quoted at rupture have no thermodynamic meaning, because any membrane under tension is inherently unstable and will be observed to rupture if the time scale of observation is sufficiently long. The experimental values quoted apply over some, unspecified, laboratory time scale. On this point we add that, as in experiment, we found many of our bilayers to rupture over the time we observed them, but the time scale for this to happen was significantly greater than that for fusion. If the bilayer area expansion, or equivalently, its tension were reduced, either in experiment or in our simulation, the time scales for the bilayers to fuse and later to rupture would both increase, perhaps to the extent of making impossible the observation of fusion. Indeed we chose the value of tension in the simulation such that fusion could be observed conveniently. One could still ask whether, in addition to increasing the time scale for fusion, a significant reduction in bilayer tension would favor an alternative fusion pathway. To attempt to answer this question, one could contemplate even longer Monte Carlo runs on membranes under less tension.

There are other physical parameters which might affect the fusion pathway but which are not encompassed by the quantities in Table I. For example, one might ask whether the fusion pathway is expected to be the same for large virus-encapsulating endosomes as it is for small synaptic vesicles. Thus one would consider the dimensionless ratio of the membrane’s hydrophobic thickness to the radius of the vesicle in question. We have considered the simplest case of planar membranes for which this ratio is zero. For endosomes encapsulating influenza viruses with an average diameter of $`100nm`$, the ratio is small, less than, but of the order of, 0.03, but for synaptic vesicles of typical diameter $`50nm`$, it is at least twice this. It is not difficult to imagine that for a sufficiently large value of this ratio, which implies a small area of contact between the fusing vesicles, there might be insufficient room for the growth and movement of the stalk we have observed, so that our mechanism would be supplanted by another. But we do not know this.

Ultimately the most meaningful test of the applicability of our mechanism to biological fusion is comparison to experiment, and our scenario does have testable consequences. First, because of the initial stalk formation, one expects to see the mixing of lipids in the two proximal layers before the fusion pore opens, if it forms at all, a result which is in accord with experiment Melikyan et al., (1995); Lee and Lentz, (1997); Evans and Lentz, (2002). Second, due to the formation of holes in each bilayer near a stalk, our scenario allows for the mixing of those lipids in the cis and trans leaves of one bilayer and also of lipids in the cis leaf of one bilayer with those in the trans leaf of the other. The standard hemifusion mechanism does not permit either process. Note that this movement is different from lipid flip-flop which is known to be very slow. Mixing of lipids between the cis and trans monolayers has been observed in fusion Lentz et al., (1997); Evans and Lentz, (2002), but it has not yet been determined from which membrane they originate and in which membrane they terminate. We have monitored the amphiphiles to see whether they remain in the leaf in which they were situated at the beginning of the Monte Carlo run, or mix with amphiphiles in other leaves. Instantaneous assignment of amphiphiles to a respective monolayer was determined by the center of mass of their hydrophilic part. The results are shown in Fig. 10. They share with experiment the fact that the membrane of origin is not distinguished nor is the membrane of final residence. In order to evaluate the results for the apposed bilayers under tension, we have also included those for the single isolated bilayer under zero tension and under the same tension ($`\gamma /\gamma _0=0.75`$) as in the simulations of fusion. Lateral tension greatly enhances the flip-flop rate in the single bilayer system. This effect can be explained by an overall thinning of the membrane, which lowers the translocation barrier, as well as by the diffusion of amphiphiles through the transient holes formed under tension. In the simulations of the apposed bilayers, translocation of amphiphiles from the trans leaves initially follows the same dynamics as in the single bilayer system, but eventually deviates from it, apparently due to the formation of holes facilitated by the appearance of the stalks, as discussed in the previous section. Amphiphiles from the cis leaves undergo mixing to the largest extent, as would be expected due to stalk formation. Third, our mechanism allows for transient leakage during fusion. As noted earlier, there will be greater leakage if fusion occurs via pathway two in which the stalk aligns and surrounds two holes than if it occurs via pathway one in which the stalk rapidly surrounds one hole before the second appears. Clearly the amount of leakage depends on the size of the transient holes formed in the bilayer, the time between the formation of the initial stalk and the completion of the fusion pore, and the diffusion constant of the molecules which leak. This constant introduces another time scale whose magnitude, relative to that of fusion pore formation, determines whether the fusion process is observed to be leaky or tight.

It is clear that within our mechanism, leakage via transient holes and fusion via pore formation are correlated in space and time. The latter is shown in Fig. 11 which presents, as a function of time, the area of holes and that of fusion pores from one of the simulation runs. One sees in this figure, as in the Monte Carlo snapshots, that the rate at which holes appear, and therefore the rate at which leakage should occur, increases significantly before, and is correlated with, the formation of fusion pores. . Once the fusion pore has formed, the creation of other holes decreases due to release of tension initially stored in the membranes.

The question of whether transient leakage is characteristic of membrane fusion is an open one. On the one hand, some experiments detect no leakage Spruce et al., (1991); Tse et al., (1993); Smit et al., (2002), while on the other there is a great deal of evidence that fusion of biological membranes is, indeed, a leaky process Shangguan et al., (1996); Dunina-Barkovskaya et al., (2000); Bonnafous and Stegmann, (2000); Haque and Lentz, (2002); Smit et al., (2002). It could be argued that observed leakage is due to the presence, in these experiments, of fusion proteins, such as hemagglutinin, which are certainly present in the vicinity of fusion, and which are known to undergo conformational changes in which part of the protein inserts itself into the target vesicle. In support of this view, one could cite the well known ability of fusion peptides to initiate erythrocyte hemolysisNiles et al., (1990). Such peptides are not included in our model. This argument is vitiated, however, by the observation that leakage is also detected in the fusion of model membranes without such peptides Lentz et al., (1997); Cevc and Richardsen, (1999); Evans and Lentz, (2002). In these experiments, large molecules such as polyethyleneglycol, are used to bring the fusing vesicles together. It would be difficult to argue that these molecules, which undergo no conformational change, are responsible for the leakage as they generate an attractive osmotic force between the vesicles precisely because their large size makes it difficult for them to enter the region where the vesicles are closely opposed.

One test that might distinguish whether leakage simply accompanies fusion or is causally related to it is provided by the observation above that in our mechanism transient leakage is correlated in time and space with fusion. Just such an experiment to measure these correlations has been carried out recently Frolov et al., (2002), and is reported in the paper accompanying this manuscript. They observe that leakage is, in fact, correlated spatially and temporally with the process of fusion. Indeed, their results comparing the time sequence of the electrical conductance arising from leakage with that arising from fusion, shown in their Fig. 5 displays a remarkable similarity to our results comparing the time sequence of the areal fraction taken up by holes with that taken up by fusion pores, our Fig. 11.

While the congruence between the predictions of our model and experiment are very encouraging, there are further tests we should like to apply to it. Foremost among these is to determine the free energy barriers for the various steps along the fusion pathway. As noted above, it is relatively simple to determine the internal energy during the course of the simulation as one need only monitor the interactions between all segments. But the simulations cannot easily evaluate the entropy changes along the fusion pathway or, therefore, the free energy barrier. To determine the actual value for the free energy barrier, calculations using self-consistent field theory, which have been extremely successful in describing the phase behavior of amphiphiles Matsen and Schick, (1994); Matsen and Bates, (1996) are currently being pursued by us. Also, elastic constants of the simulated amphiphilic monolayers, e.g. calculated in Müller and Gompper,, 2002, could be employed in the simpler phenomenological theories, which have proved to be so useful. Comparison with the full self-consistent field calculations would permit determination of the accuracy of these elastic models in describing the highly curved intermediates involved in the fusion reaction. Furthermore, there is an extensive experimental evidence on the effect of lipids of differing architecture on fusion Chernomordik, (1996); Zimmerberg and Chernomordik, (1999). The self-consistent field theory is able to describe such differences Matsen, (1995); Li and Schick, (2000) and to determine both the spatial distribution of different amphiphiles in inhomogeneous structures such as the stalk, the holes and the fusion pore, as well as the change in the free energy of these structures. Results of these investigations will be published separately.

It would be of great interest to repeat our simulations under different membrane tension, as this would help to clarify the importance of fusion peptides in bringing about such tension. Finally, it would be desirable to carry out simulations in which fusion peptides are included explicitly. One could investigate whether the membrane perturbations associated with such inclusions provide sites for the nucleation of the small holes that are necessary for the formation of the fusion pore. If this were so, one could test the further inference that, by providing nucleation sites in close proximity, one in each membrane, such peptides facilitate sucessful and rapid fusion thereby reducing leakage.

###### Acknowledgements.

We acknowledge very useful conversations with L. Chernomordik, F. Cohen, M. Kozlov, B. Lentz, D. Siegel, and J. Zimmerberg. We are particularly grateful to V. Frolov for sharing his knowledge and expertise with us. Financial support was provided by the National Science Foundation under grants DMR 9876864 and DMR 0140500 and the DFG Bi 314/17 . Computer time at the NIC Jülich, the HLR Stuttgart and the computing center in Mainz are also gratefully acknowledged. M.M. thanks the DFG for a Heisenberg stipend.


text:
# Number of loops of size ℎ in growing scale-free networks.

## Abstract

The hierarchical structure of scale-free networks has been investigated focusing on the scaling of the number $`N_h(t)`$ of loops of size $`h`$ as a function of the system size. In particular we have found the analytic expression for the scaling of $`N_h(t)`$ in the Barabási-Albert scale-free network. We have performed numerical simulations on the scaling law for $`N_h(t)`$ in the BA network and in other growing scale free networks, such as the bosonic network and the aging nodes network. We show that in the bosonic network and in the aging node network the phase transitions in the topology of the network are accompained by a change in the scaling of the number of loops with the system size.

PACS numbers: 89.75.Hc, 89.75.Da, 89.75.Fb

Natural and social systems ranging from protein-protein interactions  to the World Wide Web  can be represented as scale-free (SF) networks , that is, as a set of nodes connected by links with special statistical properties. These systems are characterized by large fluctuations in the individual degrees, i.e. in the number of links pointing or leaving each node. In a SF network, the statistical distribution $`P(k)`$ of the degree $`k`$ is a power law $`P(k)k^\gamma `$. Accordingly, if $`\gamma <3`$, the individual degree distribution displays infinite variance in the limit of infinite network size.

Besides this striking property, models traditionally defined on regular lattices display peculiar effects when defined on a scale-free topology instead . Such phenomena have increased the interest in finding universal mechanisms to generate SF networks, which seem to arise independently in so diverse contexts.

Recently, these properties have been reproduced by the Barabási-Albert (BA) dynamical model of a random network. Its simple and natural algorithm, based on a “preferential attachment” rule , triggered an avalanche of research activity in the field and generated a rich “zoology” of SF network models sharing the same fundamental ingredients. Though recently alternative mechanisms  for the generation of scale-free network have been proposed, in this work we restrict our attention to growing scale-free networks with “preferential attachment”. But the degree distribution is not the only relevant topological quantity. For the characterization of networks it is necessary to look at the motifs recurrent on it and in particular at the number of loops of length larger then three . In this work we derive a general formula for the number of loops of size $`h`$ that generalizes the known result concerning triangular loops  ($`h=3`$) for a BA network. In particular, we find that the number of loops $`N_h(t)`$ of size $`h`$ at time $`t`$ in a BA network scales as $`N_h(t)(m/2\mathrm{log}(t))^\psi `$. In addition, we show that our formula is robust when tested on other growing SF network models, such as the bosonic network model (BN)  and the aging nodes model (AN) . In the bosonic network the scaling of the number of loops with the system size changes below the Bose-Einstein phase transition becoming a power-law scaling of the type $`N_h(t)t^\xi `$.

The BA model  was the first and simplest algorithm generating scale-free undirected networks. In this model, a new node is added to the network at each time step, and it is connected by a fixed number of links $`m`$ to highly connected existing nodes (preferential attachment). According to this rule, the probability $`\mathrm{\Pi }_i`$ that an existing node $`i`$ at time $`t`$ acquires a new link is assumed to be proportional to its degree $`k_i(t)`$. This reads

$$\mathrm{\Pi }_i=m\frac{k_i(t)}{_jk_j(t)}.$$
(1)
The model can be analyzed by a mean field approximation . By this approach, one finds that the average degree of a node $`i`$ that entered the network at time $`t_i`$ increases with time as a power-law

$$k_i(t)=\frac{m}{2}\sqrt{\frac{t}{t_i}}.$$
(2)
A network built in this way displays a power law degree distribution $`P(k)k^3`$. In order to investigate the inhomogeneous topology of the network, we define a loop of size $`h`$ (a $`h`$-loop) as a closed path of $`h`$ links that visits each intermediate node only once.

In the case $`m>2`$, the BA scale-free network is a very compact network, with loops of any size. As the network evolves, new loops are introduced in the network. By definition, new loops include the newly added node: indeed, a new $`h`$-loop is formed if the new node is connected to two nodes already connected by a self-avoiding path of size $`h2`$. We indicate with $`p_{i,k}`$ the probability that the nodes $`i,k`$, attached to the network at time $`t_i,t_k`$, are connected by a link. The rate at which new loops of length $`h`$ are formed is given by the probability that the new node $`k`$ is linked to two existing nodes $`i`$ and $`j`$ times the probability $`P_{i,j}^{h2}(t)`$ that they are already connected by a self-avoiding path of size $`h2`$. Therefore, we write the following rate equation for the average number of $`h`$-loops $`N_h(t)`$

$$\frac{<N_h(t)>}{t}=\frac{1}{2}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}p_{i,k}p_{j,k}P_{ij}^{h2}(t),$$
(3)
where the factor $`\frac{1}{2}`$ takes into account that each pair of nodes $`i`$, $`j`$ has been counted twice in the sums.

On the other hand, two nodes belong to a $`h`$-loop if there is a link between them and, besides it, they are connected by a self-avoiding path of length $`h1`$. Let $`P_{i,j}^{h1}(t)`$ be the probability that this path exists; thus, the probability that the link between node $`i`$ and $`j`$ belongs to a $`h`$-loop is given by $`p_{i,j}P_{i,j}^{h1}(t)`$. We obtain the average number $`<N_h(t)>`$ of $`h`$-loops in the system times $`2h`$ by summing this quantity over all the nodes $`i`$, $`j`$ in the network. In fact, each loop has been counted $`2h`$ times, because there are $`h`$ nodes in the loop, and two possible directions. Hence, we can write

$$<N_h(t)>=\frac{1}{2h}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}p_{i,j}P_{ij}^{h1}(t).$$
(4)
By replacing Eq. $`(\text{2})`$ in Eq. (1), the probability that the node $`i`$ is attached to node $`k`$ is

$$p_{i,k}=\frac{m}{2}\frac{1}{\sqrt{t_it_k}}.$$
(5)
Moreover, the probability that a node $`k`$ is connected with nodes $`i`$ and $`j`$ is proportional to the probability that $`i`$ and $`j`$ are already connected, i.e.

$$p_{i,k}p_{j,k}=m\frac{1}{2t_k}p_{i,j}.$$
(6)
By replacing this result in $`(\text{3})`$ and by the definition $`(\text{4})`$, we obtain

$$\frac{<N_h(t)>}{t}=\frac{m}{2t}(h1)<N_{h1}(t)>.$$
(7)
Consequently, the rate at which new loops of size $`h`$ are introduced in the system is proportional to the mean number of loops of size $`h1`$.

Equation $`(\text{7})`$ has a recursive structure that allows its integration without any detailed information about the probabilities $`P_{i,j}^h(t)`$. In fact, the rate at which new loops of length $`h`$ are formed can be expressed only in terms of the number of loops of minimal size (i.e. $`h=3`$),

$$\frac{^{h3}<N_h(\zeta )>}{\zeta ^{h3}}=(h1)!<N_3(\zeta )>$$
(8)
with $`\zeta =\frac{m}{2}\mathrm{log}(t)`$. The number of triangular loops $`<N_3(\zeta )>`$ can be computed directly, for the triangular loops increase in time following $`(\text{3})`$,

$`{\displaystyle \frac{<N_3(t)>}{t}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N}{}}}p_{i,k}p_{j,k}p_{i,j}.`$ (9)
Since $`p_{i,j}`$ is given by eq. $`(\text{5})`$, approximating sums by integrals we write the rate equation in the form

$`{\displaystyle \frac{<N_3(t)>}{t}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{m}{2}}\right)^3{\displaystyle _0^t}𝑑t_i{\displaystyle _0^t}𝑑t_j{\displaystyle \frac{1}{t_i}}{\displaystyle \frac{1}{t_j}}{\displaystyle \frac{1}{t}}`$ (10)
$`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{m}{2}}\right)^3{\displaystyle \frac{1}{t}}[\mathrm{log}(t)]^2.`$ (11)
Integrating $`(\text{11})`$ we find, in agreement with ,

$$<N_3(t)>=\frac{1}{3!}\left[\frac{m}{2}\mathrm{log}(t)\right]^3.$$
(12)
Using eq. $`(\text{12})`$ in eq. $`(\text{8})`$, we compute the number of loops of size $`h`$, $`<N_h(t)>`$ and we find

$`<N_h(t)>`$ $`=`$ $`\left[{\displaystyle \frac{m}{2}}\mathrm{log}(t)\right]^h(1+O(\zeta ^1)).`$ (13)
$``$ $`\left[{\displaystyle \frac{m}{2}}\mathrm{log}(t)\right]^h.`$ (14)
The expression for the scaling of $`N_h(t)`$ with the system size $`t`$ in a BA network does not suggest a practical way to measure $`N_h(t)`$. To this purpose, one has to study the symmetrical adjacency matrix $`a`$ of the network, whose generic element $`a_{ij}`$ is defined by $`a_{ij}=1`$ if $`i`$ and $`j`$ are connected and $`a_{ij}=0`$ otherwise. Knowing this matrix, one directly measures the number of paths starting from a node $`i`$ and returning on it after $`h`$ steps that visit intermediate nodes only once. According to this argument, the term $`N_h(t)`$ has a dominating term of the type $`_i(a^h)_{i,i}/(2h)`$ and subdominant terms excluding all trivial contributions coming from paths intersecting on themselves. Let us assume that the network does not contain self-loops, i.e. $`a_{ii}=0`$ for all $`i`$ of the network. In this case, for $`h=3`$ we simply have

$$N_3=\frac{1}{6}\underset{i}{}(a^3)_{ii}.$$
(15)
For $`h=4`$ and $`h=5`$, by simple arguments it is possible to show that

$$N_4=\frac{1}{8}\left[\underset{i}{}(a^4)_{ii}2\underset{i}{}(a^2)_{ii}(a^2)_{ii}+\underset{i}{}(a^2)_{ii}\right]$$
(16)
and that

$$N_5=\frac{1}{10}\left[\underset{i}{}(a^5)_{ii}5\underset{i}{}(a^2)_{ii}(a^3)_{ii}+5\underset{i}{}(a^3)_{ii}\right].$$
(17)
Using relations $`(\text{15})`$, $`(\text{16})`$, $`(\text{17})`$, we can directly measure $`N_h(t)`$ for $`h=3,4,5`$ in the BA scale-free network model to check our analytical results.

In Fig. 1 we show that the scaling of $`N_h(t)`$ for a BA network with $`m=2`$ follows

$$N_h(t)\left(\frac{m}{2}\mathrm{log}(t)\right)^{\psi (h)}.$$
(18)
The measured effective exponent $`\psi (h)`$ reported in the Inset of Fig. 1 grows with $`h`$ as expected. Nevertheless, $`\psi (h)`$ differs from the predicted asymptotic behavior $`\psi (h)=h`$ in a significant way. This can be explained by the fact that we are considering networks of size up to $`t=10^4`$ nodes and we are still far form the asymptotic behavior $`t\mathrm{}`$.

To check the robustness of the scaling relations of $`N_h(t)`$, we have measured the number of loops of size $`h=3,4,5`$ in two alternative growing SF network models: the bosonic network (BN)  and the aging nodes network (AN) .

In the BN bosonic network, each node $`i`$ is assigned an innate quality, represented by a random ’energy’ $`ϵ_i`$ drawn form the probability distribution $`p(ϵ_i)`$. The attractiveness of each node $`i`$ is then determined jointly by its connectivity $`k_i`$ and its energy $`ϵ_i`$. In particular, the probability that node $`i`$ acquires a link at time $`t`$ is given by

$$\mathrm{\Pi }_i=\frac{e^{\beta ϵ_i}k_i(t)}{_je^{\beta ϵ_j}k_j(t)},$$
(19)
i.e. low energy, high degree nodes are more likely to acquire new links. The parameter $`\beta =1/T`$ in $`\mathrm{\Pi }_i`$ tunes the relevance of the quality with respect to the degree in the acquisition probability of new links. Indeed, for $`T\mathrm{}`$ the probability $`\mathrm{\Pi }_i`$ does not depend any more on the energy $`ϵ_i`$ and the BN model reduces to the BA model. On the other hand, in the limit $`T0`$ only the lowest energy node has non zero probability to acquire new links. Moreover, in Ref. it has been shown that the connectivity distribution in this network model can be mapped on the occupation number in a Bose gas. According to this analogy, one would expect a corresponding phase transition in the topology of the network at some temperature value $`T_c`$.

In fact, for energy distributions such that ($`p(ϵ)0`$ for $`ϵ0`$), one observes a critical temperature $`T_c`$. For $`T>T_c`$ the system is in the “fit-get-rich”(FGR) phase, where nodes with lower energy acquire links at a higher rate that higher energy nodes, while for $`T<T_c`$ a “Bose-Einstein condensate”(BE) or “winner-takes-all” phase emerges, where a single nodes grabs a finite fraction of all the links. We simulated this model assuming

$`p(ϵ)=(\theta +1)ϵ^\theta `$ and $`ϵ(0,1)`$ (20)
where $`\theta =1`$. For this distribution, it has been shown that $`T_c(\theta =1)=0.7`$ . In particular, in Fig. 2 we report $`N_h(t)`$ for $`h=3,4,5`$ as a function of the number of nodes $`t`$ in the network, at temperature $`T=1.5`$ and $`T=0.5`$ respectively.

According to simulations, in correspondence with the phase transition, the scaling of $`N_h(t)`$ drastically switches from a BA-like behavior

$`N_h(t)\left({\displaystyle \frac{m}{2}}\mathrm{log}(t)\right)^{\psi (h)}`$ for $`T>T_c`$ (21)
to a power-law scaling

$`N_h(t)t^{\xi (h)}`$ for $`T<T_c.`$ (22)
We have measured the scaling of $`N_h(t)`$ with the system size $`t`$ also for a growing network with aging nodes (AN) introduced in . The model has been motivated by the observation that in many real networks, e.g. the scientific citations network, old nodes are less cited than recent ones. In the goal of representing this feature, the probability $`\mathrm{\Pi }_i`$ to attach a link to a node $`i`$ arrived in the network at time $`t_i`$ is modified to be

$$\mathrm{\Pi }_i=\frac{(tt_i)^\alpha k_i(t)}{_j(tt_j)^\alpha k_j(t)}.$$
(23)
where $`\alpha `$ is an external parameter. As in the BA model, a new node is connected to $`m`$ existing nodes. The resulting structure of such a network strongly depends on the constant $`\alpha `$. For $`\alpha <1`$, the degree distribution is a power law $`P(k)k^\gamma `$ with an exponent monotonically increasing from $`\gamma =2`$ in the limit $`\alpha \mathrm{}`$ to $`\gamma \mathrm{}`$ in the limit $`\alpha 1`$: on the other hand, for $`\alpha >1`$ no power law is observed in the degree distribution. Therefore, this model reproduces a SF network only in the region $`\alpha >1`$.

Moreover, as observed in , in the limit $`\alpha \mathrm{}`$ the oldest node is connected to an increasing fraction of all the links, reminding the condensation observed in a bosonic network. To take into account this phenomenon, we introduce a value $`\alpha ^{}`$ such that for $`\alpha <\alpha ^{}`$ the fraction of links attached to the most connected node exceeds a finite threshold $`F`$. We expect the scaling of $`N_h(t)`$ to be different in the three regions $`\alpha >1`$, $`\alpha (\alpha ^{},1)`$ and $`\alpha <\alpha ^{}`$.

We measured the total fraction of links $`k_{max}/mt`$ attached to the oldest node in a network with $`m=2`$ and $`t=10^4`$ nodes. The threshold has been fixed at $`F=0.1`$, in order to distinguish the ’condensate’ phase from the simple scale-free phase. The value for $`\alpha ^{}`$ was found to be $`\alpha ^{}=1`$. We then measured the number of loops of size $`h=3,4,5`$ for networks made of up to $`10^4`$ nodes in the three ranges of value of $`\alpha `$. We have observed that, for $`\alpha >1`$, the number of loops of size $`h`$ scales linearly with $`t`$ (at least for $`h=3,4,5`$). In the Inset (a) of Fig. 3 we report the data for $`\alpha =1.5`$. On the contrary, for $`\alpha (\alpha ^{},1)`$ we measured the scaling

$$N_h(t)\left(\frac{m}{2}\mathrm{log}(t)\right)^{\psi (h)}$$
(24)
with $`\psi (h)`$ a monotonic function of $`h`$. In Inset (b) of Fig. 3, data for $`\alpha =0.5`$ are reported. Finally, in the region $`\alpha <\alpha ^{}=1`$, $`N_h(t)`$ becomes proportional to a power-law of the system size, as it is shown in Inset (c) of Fig. 3, referring to the case $`\alpha =5`$.

In conclusion, we have introduced the number of $`h`$-loops $`N_h(t)`$ for a network of $`t`$ nodes as a characterizing quantity for random networks. We have observed that in the BA scale-free networks $`N_h(t)`$ scales as a power law of the logarithm of the system size. Moreover, we have observed that indeed this scaling seems to be a marking feature of growing scale-free networks with preferential attachment. In particular, topological phase transitions in the bosonic network and in network including aging nodes are accompanied by a drastic change in the scaling of $`N_h(t)`$ with the system size.

The authors are grateful to Yi-Cheng Zhang and Alessandro Vespignani for their useful comments. This work has been partially supported by the European Commission - Fet Open project COSIN IST-2001-33555.


text:
# Decay rates and survival probabilities in open quantum systems

## Abstract

We provide the first statistical analysis of the decay rates of strongly driven 3D atomic Rydberg states. The distribution of the rates exhibits universal features due to Anderson localization, while universality of the time dependent decay requires particular initial conditions.

PACS numbers: 05.45.Mt, 32.80Rm, 42.50Hz, 72.15Rn

The macroscopic transport properties of classical and quantum systems sensitively depend on the dynamics of their microscopic constituents. Whereas regular Hamiltonian dynamics implies the quasiperiodic confinement of phase space density, chaotic motion and disorder generically induce diffusion-like probability transport exploring all accessible phase space at sufficiently long times. In low-dimensional, disordered quantum systems Anderson localization efficiently inhibits diffusive transport through quantum interference, as first predicted by Anderson for the charge transfer across disordered solid state samples. Later on it was realized that chaos can substitute for disorder in Anderson’s scenario, leading to dynamical localization  of quantum probability transport. However, as a further complication, quantum systems with classically non-integrable Hamiltonian dynamics (or, in the absence of a clear classical analog, with strongly coupled degrees of freedom) most often do not feature purely chaotic but rather mixed regular chaotic dynamics: besides the chaotic component of phase space the latter also comprises regular regions separated from the chaotic domain by fractal structures like cantori . Hence, dynamically localized quantum transport may be amended by quantum tunneling from regular regions or by semiclassical localization in the vicinity of partial phase space barriers .

A rather robust measure of transport is the transmission or decay probability of some initial probability distribution across or from a confined region $`\mathrm{\Omega }`$ of phase space. Whereas the transmission problem immediately suggests a scattering approach, decay is most conveniently described by the norm $`𝒫_\mathrm{\Omega }U|\varphi _0`$ of the time evolved initial state $`|\varphi _0`$ upon projection $`𝒫_\mathrm{\Omega }`$ on the phase space domain we want to characterize. If $`\mathrm{\Omega }`$ is confined by an absorbing boundary, attached to some leads which allow for transport to infinity, or coupled to some continuum of states, then $`𝒫_\mathrm{\Omega }U|\varphi _0`$ will decrease monotonically in time, for generic $`|\varphi _0`$ . Given the spectral decomposition of the Hamiltonian, this finite leakage from $`\mathrm{\Omega }`$ can be accounted for by associating non-vanishing (exponential) decay rates $`\mathrm{\Gamma }_j`$ with eigenstates $`|ϵ_j`$ with non-vanishing support on $`\mathrm{\Omega }`$. Given some initial $`|\varphi _0=𝒫_\mathrm{\Omega }|\varphi _0`$, the survival probability within the domain $`\mathrm{\Omega }`$ boils down to

$$P_{\mathrm{surv}}(t)=𝒫_\mathrm{\Omega }U|\varphi _0=\underset{j}{}w_j\mathrm{exp}(\mathrm{\Gamma }_jt).$$
(1)
All information on the decay process is now encoded  (i) in the set of decay rates – which is a global spectral property of the problem – and (ii) in the expansion coefficients $`w_j=|\varphi _0|ϵ|^2`$ of the initial wave packet in the eigenbasis, what is a specific representation of the initial condition, and a local property to the extent that $`|\varphi _0`$ is localized. With these premises, our initial assertion on the crucial importance of microscopic dynamics for macroscopic transport leads to the guiding question of our present contribution : How do classically mixed regular-dynamics affect these global and local spectral properties, and, consequently, the time dependence of $`P_{\mathrm{surv}}`$? Whilst, recently, a considerable corpus of literature has addressed this problem within the context of simple models of mesoscopic transport , we shall here underpin our rather general answer by a numerically exact treatment of highly excited three dimensional Rydberg states under strong microwave driving . These are paradigmatic objects to (theoretically and experimentally) study  quantum probability decay in the presence of tunneling as well as semiclassical and dynamical localization, and allow a natural and clean definition of $`P_{\mathrm{surv}}`$ and the associated projection $`𝒫_\mathrm{\Omega }`$, without any approximations.

Let us start with a brief reminder of dynamical localization in periodically driven hydrogen Rydberg atoms, to establish the analogy between the general transport problem depicted above and our specific choice of its physical realization: The atoms are initially prepared in an unperturbed atomic eigenstate $`|\varphi _0=|n_0\mathrm{}_0m_0`$ labeled by the principal quantum number $`n_0`$, the angular momentum $`\mathrm{}_0`$, and the angular momentum projection $`m_0`$. For linearly polarized coherent driving (which we treat in dipole approximation) of amplitude $`F`$ at frequency $`\omega `$, $`m_0`$ is a good quantum number inherited from the rotational symmetry of the problem, and the atomic excitation process proceeds through the subsequent absorption and/or emission of photons from/into the field, along a sequence of near resonantly coupled states $`|n\mathrm{}`$. Provided $`F`$ is large enough compared to the attractive force exerted on the electron by the nucleus, the classical dynamics of the driven Rydberg electron initially launched along a Kepler orbit turns chaotic, leading to an efficient energy absorption also on the quantum level, which finally leads to the ionization of the atom. Hence, we define the electronic survival probability as the remaining bound state population after an interaction time $`t`$, with $`𝒫_\mathrm{\Omega }`$ the projector on the bound state component of Hilbert space. The $`\mathrm{\Gamma }_j`$ in Eq. (1) then become the ionization rates of the Floquet eigenstates $`|ϵ_j`$ of the atom in the field, and the $`\omega _j`$ the associated overlaps (averaged over one field cycle ) of $`|\varphi _0`$ with $`|ϵ_j`$. With some theoretical and numerical effort, all these quantities can be extracted from the fundamental Floquet eigenvalue problem, for typical parameter values of state of the art laboratory experiments .

The atomic ionization process is mapped on the Anderson model through the atomic localization parameter $`=\xi /\eta `$ , with $`\xi `$ the localization length in units of the photon energy. $`\eta `$ is the ionization potential of $`|\varphi _0`$, measured in units of $`\omega `$, and is completely determined by $`n_0`$, due to the angular momentum degeneracy of the hydrogen spectrum. Provided that $`n_0`$ – through its correspondence with the principal action of the classical dynamics – is chosen within the chaotic component of phase space (possibly garnished by nonlinear resonance islands and remnants of classically regular motion immerged in the “chaotic sea”), $`\eta `$ can be conceived as a measure of the extension of the domain of complex transport along the energy axis. If we identify $`\eta `$ with the sample length in the original Anderson problem and average over different realizations of a fixed value of $`1`$, quantum suppression of transport through exponentially localized eigenfunctions over the quasi-resonantly coupled states on the energy axis  implies a distribution $`P(\mathrm{\Gamma })1/\mathrm{\Gamma }`$ of the ionization rates entering Eq. (1), in accord with the current theory  on open Anderson insulators.

Statistically independent realizations of a fixed localization parameter can be generated by simultaneous variation of $`F`$ and $`\omega `$, due to the relation $`6.66F^2n_0^2/\omega ^{7/3}`$ from the original theory on dynamical localization  in periodically driven, one dimensional hydrogen atoms. Whilst this theory has no quantitative predictive power , it is nonetheless known to give qualitatively correct predictions, even for 3D Rydberg states with not too large values of the angular momentum quantum number $`\mathrm{}`$. In the sequel we therefore use this quantity such as to guide our choice of $`F`$ and $`\omega `$ values.

Fig. 1 shows the distribution of ionization rates $`\mathrm{\Gamma }`$ of real 3D hydrogen Rydberg states and of a one dimensional model atom (where the Rydberg electron is confined to 1D configuration space $`z>0`$ defined by the polarization axis, with the Coulomb singularity at the origin ) exposed to a microwave field, within a Floquet zone (of width $`\omega `$ on the energy axis ) centered around the $`n_0=70`$ (3D) and $`n_0=100`$ (1D) Rydberg manifold. The different values of $`=0.25\mathrm{}2.0`$ are realized by sampling the spectra over frequency ranges $`\omega /2\pi =35.5\mathrm{}36.1\mathrm{GHz}`$ (3D), and $`13.16\mathrm{}16.45\mathrm{GHz}`$ (1D), respectively, and adjusting $`F`$ correspondingly. Due to the dramatically enhanced spectral density of the 3D as compared to the 1D problem (a consequence of the additional degree of freedom labeled by $`\mathrm{}`$) only ten equidistant $`\omega `$-values are needed to generate the 3D data (with a total sample size of approx. 25000 eigenstates), in contrast to 500 in the 1D case (up to 100000 eigenvalues).

For small localization parameter $`=0.25`$, both, the 1D model atom as well as the real object do exhibit ionization rates distributed according to the $`\mathrm{\Gamma }^1`$ law, over approx. seven orders of magnitude from $`\mathrm{\Gamma }10^{15}`$ to $`\mathrm{\Gamma }10^9\mathrm{a}.\mathrm{u}.`$, in nice agreement with the prediction for decay from a disordered solid . However, as the localization parameter is increased (by systematically increasing $`F`$ over the entire frequency range indicated above) to $`=1.0`$, we observe a depletion of the distribution at small rates, balanced by a decrease of the slope of $`P(\mathrm{\Gamma })`$ in an intermediate range $`\mathrm{\Gamma }10^{10}\mathrm{}10^8\mathrm{a}.\mathrm{u}.`$ Such behavior is clearly incompatible with the simple picture of exponentially localized probability densities along the energy axis (due to dynamical localization of the eigenstates over the chaotic phase space component) , and suggests the presence of alternative transport mechanisms which compete with dynamically localized diffusion . As already anticipated in the introduction, such observation does not come as a surprise, since periodically driven Rydberg states exhibit a mixed regular chaotic phase space. Indeed, it is known from experimental as well as exact theoretical/numerical studies  of strongly perturbed atomic Rydberg systems, and, more recently, from mesoscopic systems with decay , that the eigenstates of such quantum systems can be classified according to their localization  (i) on regular and/or elliptic regions, (ii) along remnants of regular motion immerged in the chaotic sea (“separatrix/hierarchical states” living on “can-tori”), and (iii) in the chaotic domain of phase space. Furthermore, elliptic regions as well as remnants of invariant tori are rather robust under changes of some control parameter such as $`F`$ in our present case .

Accordingly, the localization properties as well as the decay rates of the associated eigenstates tend to be less sensitive when changing $`F`$, whereas those states localized in the chaotic domain, with their decay rates essentially determined by $``$, exhibit rapidly increasing rates as $`F`$ is increased. Hence, we attribute the deviation from the $`\mathrm{\Gamma }^1`$ law at large values of $`F`$ in Fig. 1 to the rapid increase of the decay rates of class (iii) states with $``$ (corresponding to the effective destruction of dynamical localization), whereas class (i) and (ii) states accumulate at intermediate $`\mathrm{\Gamma }`$-values at the largest $``$-values considered here (we verified this assertion by inspection of the Husimi phase space representations of the Floquet eigenstates of the 1D problem ). At small values of $``$, the latter exhibit extremely small ionization rates below $`\mathrm{\Gamma }10^{13}\mathrm{a}.\mathrm{u}.`$ (or even below $`\mathrm{\Gamma }10^{15}\mathrm{a}.\mathrm{u}.`$, this is below our numerical precision, hence the plot does not extend into that domain), what leads to isolated, extremely narrow resonances in experimental or numerical photoionization cross sections . Therefore, in this parameter regime, $`P(\mathrm{\Gamma })`$ is dominated by the decay rates of the dynamically localized class (iii) states, in agreement with the Anderson scenario. Let us also comment on the apparently distinct weight of class (i) and (ii) states in the 1D and the 3D situation, which is born out in Fig. 1: the depletion of the low-$`\mathrm{\Gamma }`$ part is less pronounced for the real atom, and equally so the region of decreased slope at intermediate $`\mathrm{\Gamma }`$-values. This is due to the contribution of high $`\mathrm{}`$ states (at fixed $`m_0`$) to the rate distribution – unavailable in 1D – which exhibit vanishing decay rates and do not contribute to the plotted distribution at small $``$. As $``$ is increased, these states contribute to $`P(\mathrm{\Gamma })`$ in the range $`\mathrm{\Gamma }<10^{11}\mathrm{a}.\mathrm{u}.`$, such as to compensate for the re-shuffling of low-$`\mathrm{}`$ states from small to intermediate $`\mathrm{\Gamma }`$’s. This smooths the 3D distribution as compared to the 1D result.

Hence, we have seen how chaotic transport can mimic signatures of disorder in the quantum decay rate distribution, and induce amendments thereof due to peculiar structures of mixed regular chaotic dynamics. As obvious from Eq. (1), $`P(\mathrm{\Gamma })`$ does not determine the (experimentally often very easily accessible ) survival probability alone, but only in conjunction with the expansion coefficients $`w_j`$. Starting from a general form $`P(\mathrm{\Gamma })\mathrm{\Gamma }^\alpha `$, many authors infer a “universal” time dependence $`P_{\mathrm{surv}}1/t^{2\alpha }`$, leaning on the additional assumption $`w_j\mathrm{\Gamma }_j`$ . The latter is nonetheless only justified if the initial state $`|\varphi _0`$, which represents a local probe of the spectrum, is localized close to the boundary of $`\mathrm{\Omega }`$ – which is the “edge of the sample” in a language adapted to the Anderson scenario, and the ionization threshold $`n\mathrm{}`$ in the atomic problem. For a generic choice of $`|\varphi _0`$ “within the sample”, i.e. somewhere within $`\mathrm{\Omega }`$, this basic assumption is utterly misleading, as illustrated in Fig. 2. In Fig. 2 a), we show the distribution of the $`w_j`$ and of the corresponding $`\mathrm{\Gamma }_j`$, on a double logarithmic plot, for $`|\varphi _0=|n_0=70\mathrm{}_0=0`$ in driven 3D hydrogen. The essentially uncorrelated $`w_j`$ and $`\mathrm{\Gamma }_j`$ nicely correlate with the earlier finding that phase space localization properties and decay rates of strongly driven Rydberg atoms are not unambiguously related , and clearly devalidate the above proportionality assumption. On the other hand, it can easily be shown within a simple, one dimensional Anderson model how the generic situation of Fig. 2 a) continuously evolves into the situation of strongly correlated $`w_j`$ and $`\mathrm{\Gamma }_j`$, as illustrated in Fig. 2 b): There, we calculated the decay rates of the eigenstates of a 1D tight binding Hamiltonian of length $`N=1000`$, with diagonal disorder and constant off-diagonal terms. Decay is introduced in this model by absorbing boundary conditions at one end, whereas a reflecting wall is imposed at the other. If $`|\varphi _0`$ is placed well inside the sample, expansion coefficients and decay rates are uncorrelated, though the distribution collapses onto a straight line $`w_j\mathrm{\Gamma }_j`$ as $`|\varphi _0`$ is shifted towards the open end (site $`n_\mathrm{i}=1000`$) of the sample. It is this latter situation which was realized in recent numerical calculations on periodically driven hydrogen atoms additionally subjected to an intense static electric field .

Given the results illustrated in Fig. 2, we can no more hope for a universal decay law for the survival probability $`P_{\mathrm{surv}}(t)`$ – despite the rather universal features of $`P(\mathrm{\Gamma })`$. This is finally illustrated in Fig. 3, where we show $`P_{\mathrm{surv}}(t)`$ for 1D and 3D hydrogen in the case of large $``$: Whilst algebraic decay of $`P_{\mathrm{surv}}(t)t^\gamma `$ is generic (even if modulated by local fluctuations ), simply due to the broad distribution of the $`w_j`$ and $`\mathrm{\Gamma }_j`$ as illustrated in Figs. 1 and 2, the “universal” decay exponent $`\gamma =2\alpha `$ is only observed in accidental cases – or if we reshuffle the expansion coefficients of Fig. 2 a) according to $`w_j\mathrm{\Gamma }_j/\mathrm{\Gamma }`$, with $`\mathrm{\Gamma }`$ the average rate over the entire spectrum, as illustrated in Fig. 3 b). In this, however, unrealistic case, $`P_{\mathrm{surv}}(t)`$ asymptotically decays with $`\gamma 1.1`$, which perfectly matches the prediction $`\gamma =2\alpha `$, with $`\alpha =0.9`$ extracted from the small $`\mathrm{\Gamma }`$ ($`<10^{10}\mathrm{a}.\mathrm{u}.`$) part of the corresponding data set in Fig. 1 b). Therefore, universal features of the statistics of the decay rates of classically chaotic quantum systems do generically not carry over to the time dependent decay of the survival probability of some initial state prepared at an arbitrary location in phase space. This simply expresses the essential decorrelation of the phase space localization properties of the eigenstates of chaotic quantum systems and of their asymptotics – which alone determine their decay.

Support from EU Program QTRANS RTN1-1999-08400, is gratefully acknowledged (S.W.).