text:
# Real event detection and the treatment of congestive heart failure: an efficient technique to help cardiologists to make crucial decisions

## Abstract

Using a method of entropic analysis of time series we establish the correlation between heartbeat long-range memory and mortality risk in patients with congestive heart failure.

This paper deals with an advanced aspect of statistical mechanics whose recent demonstration  is proven here to afford an efficient criterion to deal with the mortality risk of patients with Congestive Heart Failure (CHF). We define the condition corresponding to the mean position of healthy subjects in a plane, called, as we shall see, physiological plane. We show that the distance of a CHF patient from this optimal condition correlates with mortality, and that the CHF patients very close to it survived. The importance of this criterion does not need to be defended. We think that the demonstration of this important property is of interest for both researchers at the frontier of statistical mechanics and cardiologists. The first step of our approach is based on observing a sequence of numbers $`\{T_i\}`$, denoting the distance between two nearest neighbor pulses of a given electrocardiogram (ECG). One way to denote the time between two pulses is to measure the time elapsed betwen two adjacent R waves in the recorded electrical signal. This time is referred to as the RR interval, and the resulting sequence is consequently denoted as RR sequence.

We considered for this study 13 male CHF patients, from a study base of 320 subjects, who experienced cardiac death during a follow-up of 26 months (average 19 months, median 22 months). Inclusion criteria were absence of pulmonary or neurological disease, absence of acute myocardial infarction or cardiac surgery within the previous 6 months, absence of any other disease limiting survival, stable therapy for at least 2 weeks and good quality 24-hour Holter recordings, with an ectopy rate less than 5%. A comparable number of control subjects (16 patients), matching for age, sex, NYHA class (a functional and therapeutic classification for prescription of physical activity for cardiac patients) and etiology, was then selected. These latter patients did not experience cardiac death after follow-up. All patients had a 24-hour Holter recording at baseline, together with standard functional evaluation including measurement of left ventricular ejection fraction (LVEF), peak VO2 and Sodium (Na). Finally, RR series for 10 healthy subjects were taken from the NOnLinear TIme Series AnaLysIS (NOLTISALIS) archive. This latter data set is the result of the collaboration of several interdisciplinary Italian research centers. Experienced analysts edited these Holter recordings, manually correcting interbeat times due to ectopic beats. This editing work yielded the RR sequence, from which we generated the time series $`\{T_i\}`$.

The meaning of a given value $`T_i`$, as earlier stated, is the time distance between the $`ith`$ and the $`(i+1)th`$ pulse. The sequence $`\{T_i\}`$ can be studied as a new time-series, with $`i`$ playing the role of “time”. Moreover, the value $`T_i`$, expressed as a function of $`i`$ with $`i1`$, can be thought of as a function $`T(t)`$, namely, as a function of a continuous time variable $`t`$.

The real curve $`T(t)`$ looks erratic and disordered. However, our method of analysis shows that it is quite different from a random process. For the reader to get an intuitive understanding of this attractive but perplexing conclusion, let us describe first an ideal model, with extended memory (EM), for the time evolution of $`T(t)`$. First of all, we assume that for a given time $`\tau _{em}`$, the curve $`T(t)`$ keeps a given slope $`\alpha `$, then it abruptly gets a new slope, $`\alpha ^{}`$, for an interval of time $`\tau _{em}^{}`$, after which a new abrupt transition to a new slope $`\alpha ^{\prime \prime }`$ takes place, for a time $`\tau _{em}^{\prime \prime }`$, and so on. It is evident that the resulting $`T(t)`$ has the form of a zig-zag curve. We shall refer to the individual straight line intervals of this curve as *laminar regions*. Any laminar region is associated with its own $`\tau _{em}`$. Then we introduce the extended memory property. This is done by assuming for the waiting time distribution $`\psi (\tau _{em})`$ the following inverse power law form

$$\psi (\tau _{em})=(\nu 1)\frac{\left[\tau _{em}(\nu 1)\right]^{\nu 1}}{\left[\tau _{em}(\nu 1)+\tau _{em}\right]^\nu },$$
(1)
with $`2<\nu <3`$, where $`\tau _{em}`$ is the average waiting time. This means that if the EM model is directly observable, we can derive from $`T(t)`$ the sequence $`\{\tau _{em}(j)\}`$, where the discrete index $`j`$ denotes the time order of a given laminar region.

It is interesting to notice that this dynamic process is essentially equivalent to the strong anomalous diffusion model recently proposed by the authors of Refs. to explain the effects of a ballistic mechanism in the Rayleigh-Bénard convection. This model, in turn, is nothing but a generalization of the dynamic approach to Lévy statistics proposed years ago by the authors of Ref. . In fact, the model of Ref.  is recovered from the model we are adopting here, by assuming that the slope $`\alpha `$ has only two values, of equal intensity and opposite sign . With this equivalence in mind, we adopt the specific walking prescription proposed by Ref.. We consider a given time $`t`$ and we evaluate the number of laminar regions that have been completed within this time. Let us call this number $`n`$. Then the trajectory is

$$y(t)=nW.$$
(2)
This means that the random walker makes a jump ahead, by the same quantity $`W`$, at the end of any laminar region. The quantity $`W`$ is arbitrary and in the following we assume $`W=1`$. Then, according again to the prescriptions of Ref. , we create the trajectories

$$x(l)=y(t+l)y(t).$$
(3)
It is evident that we can move the index $`t`$ from $`t=0`$ to $`Ll`$, where $`L`$ is the total time length of the sequence under study. This makes it possible for us to evaluate the probability distribution $`p(x,l)`$, which at $`l=0`$ is a delta of Dirac located at $`x=0`$, broadening upon time increase. Note that an important step of our approach rests on the determination of the Shannon entropy of this distribution, namely

$$S(l)=_{\mathrm{}}^{\mathrm{}}𝑑xp(x,l)\mathrm{log}p(x,l).$$
(4)
This is the reason why this technique of analysis has been termed Diffusion Entropy (DE) method . According to the theory of Ref. we immediately reach the conclusion that $`p(x,l)`$ fulfils the scaling property

$$p(x,l)=(1/l^\delta )F(x/l^\delta ),$$
(5)
with $`\delta `$ being the scaling parameter, which is related to $`\nu `$ by

$$\delta =1/(\nu 1).$$
(6)
In this specific case, it is straightforward to prove, by plugging $`p(x,l)`$ of Eq. (5) into Eq.(4), that

$$S(l)=A+\delta \mathrm{log}(l).$$
(7)
This means that the DE should yield for $`S(t)`$, expressed in a linear-log plot, a straight line, whose slope is the searched value of the scaling parameter $`\delta `$.

This way of proceeding is impossible in practice, because the real $`T(t)`$ curve significantly departs from the zig-zag form of the EM model. We make the conjecture that the departure from this ideal condition is caused by the fact that the actual signal $`T(t)`$ is the superposition of the EM model signal and a much stronger, but totally random component. This makes it impossible for us to directly evaluate $`\psi (\tau _{em})`$. We proceed as follows. Let us represent the $`T(t)`$ time evolution in a $`(T,t)`$ plane, with the ordinate referring to $`T`$ and the abscissa to $`t`$. The ordinate axis is divided into cells of equal size, called $`s`$. This means that we divide the $`(T,t)`$ plane into strips of size $`s`$, and that in the ideal case of constant frequency the trajectory $`T(t)`$ would move forever remaining within the same strip. Actually, transitions from one strip to the other occur frequently. We call these transitions markers. These markers might have quite different origins. The majority of the markers are determined by the short-time noise. Many other markers correspond to the hidden laminar regions of the underlying EM model; we call these markers “pseudoevents”. As in Ref. , we indicate with the term pseudoevent a marker that does not correspond to an unpredictable transition, but it is a consequence of the division of the $`(T,t)`$ plans into strips. Only a very small number of markers coincide with the turning points of the EM signal, or are sufficiently close to them. We call these markers *real events*.

Now we have to explain why the DE is sensitive only to the real events. The time distance between the $`j`$-th and the $`(j+1)`$-th marker, defines the time $`\tau _{exp}(i)`$ of the experimental sequence $`\{\tau _{exp}(i)\}`$. It is evident that even in the case where $`T(t)`$ were an exact realization of the EM model, the waiting time distribution $`\psi (\tau )`$ might turn out to be totally different from $`\psi (\tau _{em})`$. In fact, if $`s`$ is very small the same laminar region is divided into many smaller time intervals, with the same length. These are pseudoevents. It is evident that these pseudovents do not contribute to the spreading of the distribution $`p(x,l)`$, and consequently do not contribute to the entropy increase. This means that the DE method is insensitive to pseudoevents.

What about the short-time random events? They, in principle, contribute the entropy increase, and consequently could affect the determination of the crucial parameter $`\delta `$. We can show that they do not. Notice that in the presence of the short-time random component the actual signal $`\tau _{exp}(t)`$ is given by

$$\tau _{exp}(t)=a\tau _{st}(t)+bR(t),$$
(8)
with $`a1`$ and $`b`$ close to $`1`$. The first contribution corresponds to the EM model, and the second is generated by short-time random fluctuations. The correlation function of this signal is

$$C_{exp}(t)=pC_{st}(t)+(1p)C_{random}(t),$$
(9)
where $`C_{st}(t)`$ is an inverse power law relaxation and $`C_{random}(t)`$ a relaxation function decaying to zero in one time step. If we take, without loss of generality $`\tau _{st}^{}{}_{}{}^{2}=R^2=\tau _{exp}^{}{}_{}{}^{2}`$, implying $`a^2+b^2=1`$, then we have $`p=a^2`$ and $`(1p)=b^2`$. How should parameter $`p`$ be evaluated? This is easily done monitoring the experimental correlation function at the first time step. According to the fact that $`C_{random}(t)`$ decays to zero in one step, while $`C_{st}(t)`$ is much slower, we immediately obtain $`p=C_{exp}(1)`$. Its evaluation is not independent of that of the other parameter, $`s`$. This is so because the experimental evaluation of $`\tau _{exp}(t)`$ is dependent on $`s`$.

However, while $`p`$ is $`s`$-dependent, the parameter $`\delta `$ is not. The DE method has the surprising capability of yielding a value for $`\delta `$ that is independent of $`s`$, even in the case when a strong short-time random component is present, and not only when the ideal EM model applies. How can it be so? This is so because the EM model component yields superdiffusion, while the random component generates ordinary diffusion. In the asymptotic limit of very large values of $`l`$, the superdiffusion component, which is faster than the ordinary diffusion, becomes predominant, and the DE method detects again the correct scaling $`\delta `$ even in the case where $`p1`$.

The parameter $`p`$ is very important, since it defines the statistical weight of the EM component present in the experimental signal $`T(t)`$. However, its dependence on $`s`$ makes its use questionable. However, we see that all curves $`p(s)`$ share the same properties of getting small values for both small and large values of $`s`$, with a clearly defined maximum in between, which is referred to by us as $`\pi `$. This maximum is a property independent of $`s`$. Infact, we note that almost all the healthy patients get their maximum at 30 ms, while most of the CHF patients get theirs at 20 ms. The typical dependence of $`p`$ on $`s`$ is illustrated in Fig. 1, where we can see both the case of 5 healthy subjects with $`\pi `$ located at 30 ms, and 5 CHF patients, with $`\pi `$ located at 20 ms. We think that the parameter $`\pi `$ is a reliable measure of the EM component. Consequently, we decided to represent the physiological conditions of all patients, healthy and CHF alike, in the $`(\delta ,\pi )`$ that we call *physiological plane*. The criteria adopted to define the physiological plane make the resulting “phase-space” diagram independent of the coarse-graining parameter $`s`$, and the location of any patient in the physiological plane is an objective property independent of the coarse graining parameter $`s`$. A mere visual inspection shows that the healthy and CHF patients do not mix but in relatively small region. The overlap region is so small as to make it possible for us to claim that healthy and CHF patients belong to two distinct regions of the physiological plane .

The division of the CHF patients into two groups, dead and alive, make the result of our analysis still more remarkable. To show this important property we proceed as follows. First of all, we define the center of gravity of the healthy patient, denoted with a white square in Fig. 2. We call this point optimal condition. Then, for any CHF patient we measure the Euclidean distance from the optimal condition, thereby making it possible for us to rank the CHF patients in order, according to this distance. In other words, the first CHF patient is the one with minimum distance from the optimal condition. Then we observe the remaining patients, and we rank as second the one with minimum distance from the optimal condition, and so on. We find that the first 7 patients are alive. The eighth patient is dead, and from now on the patients are either alive or dead. This suggests that the closer the patient to the optimal condition the higher the survival probability.

To support in a more rigorous way this important property, we apply the Mann-Whitney method . This is a non-parametric test, namely it does not rest on Gaussian distributions for the data. Let us count the number of patients dead, $`N_{dead}`$, and the number of patients alive, $`N_{alive}`$. The data of Fig. 2 refer to a case where $`N_{dead}=13`$ and $`N_{alive}=16`$. Let us consider the group with the smaller number of individuals. This means the group of dead patients. Then, let us evaluate the sum of the ranks of this group and let us call it $`\mu _{exp}`$. We note that $`\mu _{exp}=246`$. Under the hypothesis of no correlation between our parameters and the death probability, this resulting sum has a probability distribution that, for more than 8 elements is expected  to be approximately a Gaussian with mean $`\mu _T`$ and standard deviation $`\sigma `$ given by

$$\mu _{dead}=N_{dead}\frac{N+1}{2}\sigma =\left(N_{dead}N_{alive}\frac{N+1}{12}\right)^{1/2}.$$
(10)
From the data of Fig. 2 we obtain $`\mu _{dead}=195`$ and $`\sigma =22.8`$, while, as earlier noticed, $`\mu _{exp}=246`$. This means that the distance of $`\mu _{dead}`$ from $`\mu _{exp}`$ is $`2.37\sigma `$. In practice, we are allowed to rule out the hypothesis of no correlation between the death of CHF patients and their distance from the optimal condition: the probability for the distribution of dead and alive patients of Fig. 2 to be fortuitous, is less than 3%.

Moreover, it is important to remark that the alive patients corresponding to points in the physiological plane far from the optimal conditions, either had a serious pathology, being classified as NYHA class III (severe physical limitations, they are confortable only at rest) and therefore required a heart transplant anyway, or had a very short follow-up time (less than 6 months). Only six alive patients did not belong to either of the above conditions, but it is remarkable that all six of them occupy positions which overlap with the zone of the healthy subjects. Unfortunately, the small number of sequences available to us at this time does not allow us to calculate the survival curves, or attempt any further conclusion.

It is convenient to compare the result of this paper to the literature in the field of nonlinear or fractal analysis of chardiological data. The most recent examples are those of Refs.. Although based on different perspectives, multifractality  and memory beyond memory as a sign of healthy physiological condition, neither of these two groups could address the ambitious step of granting physicians a reliable criterion to make crucial decisions about the CHF patients. These findings, however suggest that $`\delta `$ and $`\pi `$ indexes, and especially the distance from the optimum condition in the physiological plane, should be considered for inclusion in the candidate predictors’ list of future large-scale prospective studies for risk stratification of CHF patients.

*Acknowledgements*. We gratefully acknowledge financial support from the Army Research Office through Grant DAAD 19-02-0037. One of the authors (P. A.) acknowledges European Commission POESIA project (IAP2117/27572) for financial support. L. P. acknowledges ENEL Research Grant 3000021047 for support. We also thank S. Maugeri Foundation for supplying the data.


format: text

text:
# Electronic resonance states in metallic nanowires during the breaking process simulated with the ultimate jellium model

## I Introduction

The miniaturization of the electronics components is of great importance in the development and improvement of new devices for applications in a wide number of fields. Although the laws of Nature are the same for macroscopic and mesoscopic systems, the miniaturisation process is achieving the limit where the quantum behavior of matter starts to play an important role.

If the size of the system under consideration is only a few nanometers, the atomic character of matter emerges and it can not be considered as a continuum. The regime of quantum behavior is reached also if one of the spatial dimensions of the system is reduced down to the Fermi wavelength of the conducting electrons. Then, the confinement splits the continuous electronic band in this direction into a set of discrete energy levels. In both cases, the behavior of the system changes from what is expected from the macroscopic case. In metallic nanowires the Fermi wavelength is of the same order of magnitude as the atomic distance, and both atomic and electronic discrete character compete and/or couple, determining the properties of nanowires.

There are many experimental and theoretical works which have gone deep into the understanding of the main features of nanowires. Experimental studies have focused on the investigation of the mechanical and electronic properties such as force, atomic structures and conductance, pointing out the close relation between them. Among the experimental setups we want to emphasize the role of the scanning tunneling microscopeRubio et al. (1996); Rubio-Bollinger et al. (2001); Brandbyge et al. (1995) (STM) and the mechanically controllable break-junction (MCBJ) techniques Yanson et al. (2000, 1999, 2001). In both techniques metallic nanowires are produced by putting two protrusions in contact and then pulling them away from each other over atomic distances. In this process, a nanowire is produced which upon pulling is elongated and narrowed until it eventually breaks. These methods have allowed the study of transport properties and stability of nanowires.

The MCBJ techniques have demonstrated the existence of electronic and atomic shell structures Yanson et al. (2000, 1999, 2001), analogous to those found in atomic clusters Brack (1993); de Heer (1993). In these experiments the conductance has been studied by building histograms of the conductance during the breaking process. The results show that there are conductance values which are much more probable than others. Due to the relation between the conductance and the radius at the narrowest part of the nanowire, that means that there are magic radii with enhanced stability while other radii are less stable, and therefore they appear less frequently in the conductance histograms. The atomic structures of nanowires in the last steps before breaking have been studied also with these techniques Rodrigues et al. (2000); Rubio-Bollinger et al. (2001); Rodrigues et al. (2002); Ohnishi et al. (1998); Yanson et al. (1998).

The experiments discussed above have been accompanied by supporting theoretical investigations, which can be split in two groups. The first group includes classical and ab initio molecular dynamics simulations, in which the atomic structure of nanowires is taken into account. These investigations have been successful in many aspects, e.g., showing the atomistic mechanisms of the narrowing process (appearance of dislocations, order-disorder stages, etc.) and their link with other measurable quantities such as the elongation force or the conductance Landman (1998); Sørensen et al. (1998). Moreover, from the viewpoint of the present work we notice the predictions of special atomic arrangements in STM tips and nanowire necks Landman (1998); Häkkinen and Manninen (1998); da Silva et al. (2001); Sánchez-Portal et al. (1999); Tosatti et al. (2001); Barnett and Landman (1997). The second group of models is more related to properties due to the confinement of electrons in reduced dimensions, and ignores the atomistic structure of matter. In these calculations analytic approximations as well as self-consistent electronic-structure models have been used, mainly within the jellium framework. The results obtained with these methods are also enlightening, explaining the cohesive and electronic transport properties of nanowires, especially in the case of alkali metals with strong free-electron character Ogando et al. (2002); Puska et al. (2001); Zabala et al. (1999); Yannouleas et al. (1998); Yannouleas and Landman (1997).

The aim of this paper is to simulate the breaking of nanowires. For this purpose we choose the jellium model and the self-consistent electronic-structure calculations within the density-functional theory. In spite of their simplicity jellium models have provided a simple and transparent way to understand the physics of metallic nanowires. More specifically, we use the ultimate jellium (UJ) model. This model was first proposed by Manninen Manninen (1986) to investigate the structures of alkali metal clusters. It has been used for the same purpose also in later studies Koskinen et al. (1995a, b). To our knowledge the present work is the first time the UJ-model is used to simulate the nanowire breaking. In practice, we solve the ensuing Kohn-Sham equations in a real-space point grid using the powerful Rayleigh Quotient MultigridMandel and McCormick (1989); Heiskanen et al. (2001) (RQMG) method implemented in the program package MIKATorsti et al. (2002) (Multigrid Instead of K-spAce).

Within the UJ approach, the background positive charge density is fully relaxed in shape and density so that it equals at every point with the electron density. One can think that this freedom of the positive background charge mimics the efficient rearrangement and diffusion of ions at temperatures close to the melting point at which the shell- and supershell-structure studies by the MCBJ techniques have been performed for alkali metals Yanson et al. (1999, 2000). In principle, there is no restriction for the geometry of the constriction. This is in contrast with the previous jellium calculations, which introduced ad hoc shapes for the nanowire. In our model the electrons themselves acquire self-consistently the shape, which minimizes the Kohn-Sham energy functional, and carry along the positive background charge. However, in order to reduce computational demands and to highlight the important phenomena from the complexity of possible solutions, we restrict the shapes of nanowires to the axial symmetry, i.e. rotational invariance with respect to an axis.

One of our main results is that in the narrowest part of the nanowire, electronic cluster derived structuresBarnett and Landman (1997); Landman (1998) (CDS’s) appear. This tendency of electrons to form embedded clusters in the jellium constrictions is analogous to the preferred cluster-like arrangements of atoms in contacts described by first-principles atomistic calculations by Barnett and Landman Barnett and Landman (1997); Landman (1998). CDS’s have later been reported also by other authors Häkkinen and Manninen (1998). The main difference is that in our jellium model the atomistic character of the previous works is lost and the electrons alone are responsible for the phenomenon. The single-electron states provided by the jellium model can be studied in order to gain insight into the localization effects associated with the CDS. The conductance of the constriction can be estimated either by counting the bands crossing the Fermi level, or by using the WKB formula.

The rest of the paper is organized as follows: in Section II we describe the practical features of the UJ-model and the RQMG-method to calculate the electronic structure during the elongation process. In Section III we discuss the results for the electronic properties. As a starting point, we consider the results for infinite wires. Then we focus our attention on the breaking process of a finite cylindrical UJ nanowire supported by leads. Section IV contains the conclusions.

## II Theory

### II.1 Jellium models

The jellium model has been widely used in self-consistent electronic-structure calculations of nanostructures. It simplifies the problem by replacing the ions by a uniform rigid positive charge density background, which globally neutralizes the electron negative charge. The effective potential of the Kohn-Sham Jones and Gunnarsson (1989) equations is written as (Atomic Hartree units are used throughout this paper to write the equations)

$$V_{\mathrm{eff}}=\frac{n_{}(r^{})n_+(r^{})}{|rr^{}|}𝑑r^{}+v_{\mathrm{xc}}[n_{}(r)],$$
(1)
where the first term on the right-hand side includes the electron-electron and electron-positive background Coulomb interactions and the second term gives the exchange-correlation potential within the local density approximationPerdew and Zunger (1981); Ceperley and Alder (1980) (LDA).

Different types of jellium approaches have been introduced. The simple jellium (SJ) model has the problem that there is only one equilibrium charge density, at $`r_s4.18a_0`$ ($`n_{}=3/(4\pi r_s^3)`$) corresponding approximately to the average conduction electron density in Na or K metals. This means that for $`r_s`$ values lower (higher) than $``$4.18 $`a_0`$ the jellium system tends to expand (compress). In the SJ-model the electron density has the same mean value as the positive background due to the electrostatic forces. The SJ-model gives incorrect values for properties such as the cohesive energy, surface energy and bulk modulus, due to the trend of the system to compress or expand. To improve the results corrections can be added to the SJ-model Shore and Rose (1999), e.g. using the so-called stabilized jellium model introduced by Perdew Perdew et al. (1990) and Shore Shore and Rose (1991).

In this work we use the UJ-model, the philosophy of which differs from the stabilized jellium model in that it does not try to correct the above-mentioned deficiencies of the SJ-model. The peculiarity of the UJ-model is that the positive charge background is allowed to relax. The UJ-model represents the ultimate limit in which the positive background is completely deformed to have the same density as the electrons locally at every point. In this way, the Coulomb term in the potential always vanishes, and in Eq. (1) only the exchange-correlation term survives. The total energy is then minimized in the interplay between the exchange-correlation and the kinetic energies.

One limitation of the UJ-model is that, as in the SJ-model, there is only one equilibrium charge density, at $`r_s4.18a_0`$. But, the absence of electrostatic potential disables the mechanism to keep the electrons at a given density, and inside the UJ the mean electron density becomes equal to the equilibrium density. Another property of the UJ-model, derived also from the absence of electrostatic potential, is that the shape of the electron density is to a large extent uncontrollable, and it evolves until the ground state is achieved. This property has been used to study the most favorable shapes of simple-metal atom clustersManninen (1986); Koskinen et al. (1995a, b); Reimann et al. (1998). In the present work, however, we have to deal with open systems and we have to impose certain controlling restrictions in order to model the pulling of the nanowires. The description of the solutions to these requirements is postponed to Section III.3.

### II.2 Numerical methods

In Section III.1 infinite uniform cylindrical wires are studied. Since these systems are translationally invariant along the wire axis, the relaxation of the positive background charge and electron density is limited in the radial direction. Consequently, it is necessary to solve numerically only the radial part of the Schrödinger equation (see Zabala et al.Zabala et al. (1999) for technical details).

For the systems studied in Sections III.2 and III.3, however, the translational invariance is not required. But, in addition to the rotational invariance, periodicity in the axial direction is assumed with unit cell length $`L_{\mathrm{cell}}`$. Thus, the wave functions $`\psi `$ are indexed by the quantum numbers $`m`$, $`n`$, and $`k_z`$. Here, $`m`$ is the angular momentum quantum number and $`k_z`$ is the Bloch wave vector along the wire axis. With $`m`$ and $`k_z`$ given, $`n`$ enumerates the orthogonal states in the order of increasing energy eigenvalue. The UJ system is solved by finding the self-consistent solution to the set of equations

$$\psi _{mk_zn}(r,z,\varphi )=e^{im\varphi }U_{mk_zn}(r,z),$$
(2)
$$U_{mk_zn}(r,z+L_{\mathrm{cell}})=e^{ik_zL_{\mathrm{cell}}}U_{mk_zn}(r,z).$$
(3)
$$\frac{1}{2}\left(\frac{1}{r}\frac{}{r}+\frac{^2}{r^2}\frac{m^2}{r^2}\frac{^2}{z^2}+2V_{\mathrm{eff}}(r,z)\right)U_{mk_zn}(r,z)=\epsilon _{mk_zn}U_{mk_zn}(r,z).$$
(4)
$$n(𝐫)=2\underset{mk_zn}{}(2\delta _{0m})f_{mk_zn}|U_{mk_zn}(𝐫)|^2,$$
(5)
$$V_{\mathrm{eff}}(r,z)=V_{\mathrm{xc}}(r,z)=\frac{\delta E_{\mathrm{xc}}[n(r,z)]}{\delta n(r,z)}.$$
(6)
The effective potential $`V_{\mathrm{eff}}(r,z)`$ equals the exchange-correlation potential $`V_{\mathrm{xc}}(r,z)`$. The electron density $`n(r,z)`$ is obtained by summing single-electron densities with the occupation numbers $`f_{mk_zn}`$. The degeneracies of the states are taken into account by the factor $`2(2\delta _{0m})`$ and the occupation numbers $`f_{mk_zn}`$ obey the Fermi-Dirac statistics with the Fermi level ($`E_F`$) so that the system is neutral. A finite temperature of 1200 K is used to stabilize the solution of the set of equations.

The Schrödinger equation (4) is discretized on a regular two-dimensional $`(r,z)`$ point mesh. We use standard fourth-order central-difference discretizations for the first and the second derivatives. The grid is surrounded by a frame with the thickness of two grid points. These ghost points are necessary for the evaluation of the derivatives near the edges of the computation volume. The wave functions are required to vanish at the ghost points corresponding to the radial surface of the cylindrical computation volume, whereas at the axis the values at ghost points can be evaluated by noting that $`U(r,z)=(1)^mU(r,z)`$. The Bloch-condition (Eq. (3)) gives the recipe for obtaining the values at the ghost points of the periodic boundary.

The problem with standard real-space relaxation methods for Eq. (4) is the so called critical slowing-down phenomenon resulting from the fact that at a time they use information from a rather localized region of space. As a result of the locality, the high-frequency error, corresponding to the length scale of the grid spacing, is reduced very rapidly in the relaxation. However, once the high frequency error has been effectively removed, the very slow convergence of the low-frequency components dominates the overall error reduction rate, i.e. critical slowing-down occurs. Multigrid methods avoid this problem by treating the low frequency components of the error on coarser grids, where their wavelength is comparable to the grid spacing.

Applying the multigrid methods to the Schrödinger equation is a fairly complicated task because one has to solve both the eigenvalue and the wave function simultaneously – this makes the problem nonlinear. Also, one has to solve several wave functions simultaneously, avoiding the bottleneck of orthogonalizations as well as possible. Standard methods based on the full-approximation-storageA. Brandt and Ruge (1983) method require, that the wave functions are well representable on the coarsest grid used, implying severe limitations on the acceleration obtained by the multigrid idea. We use the recently developed generalization of the Rayleigh-quotient multigrid (RQMG) methodMandel and McCormick (1989); Heiskanen et al. (2001) as implemented in the MIKA-packageTorsti et al. (2002), that avoids the problems described above. In short, one applies the Gauss-Seidel method on the finest grid. On the coarser grids one applies coordinate relaxations on the functional

$$\frac{\psi _n|H|\psi _n}{\psi _n|\psi _n}+\underset{i=1}{\overset{n1}{}}q_i\frac{\left|\psi _i|\psi _n\right|^2}{\psi _i|\psi _i\psi _n|\psi _n}.$$
(7)
This functional, which is actually defined on the finest grid, is the sum of the Rayleigh quotient and a penalty functional, which is introduced to insure the orthogonality. Moreover, the relaxations are performed simultaneously for all wave functions. See Heiskanen et al.Heiskanen et al. (2001) for a more thorough discussion of technical details.

The Kohn-Sham equations have to be solved self-consistently. In other words one has to iterate until the output potential $`V_{\mathrm{eff}}`$ obtained from Eq. (6) equals the input potential $`V_{\mathrm{eff}}`$ that is used in Eq. (4). In typical cases of electronic structure calculations, to avoid divergence due to charge sloshing, one uses sophisticated strategies to construct the input potential for the next iteration as an optimized mixture of input and output potentials of previous iterationsD. Raczkowski and Wang (2001); Annett (1995). In the UJ-iterations, however, the output potential can be taken directly as the input potential of the next iteration resulting in a rapid convergence. This is because of the absense of the long-range Coulomb interaction, which is the cause of the charge sloshing phenomenon.

## III Results and discussion

### III.1 Infinite uniform cylindrical wires

The main results of this paper concerning the nanowire breaking process are discussed in Sections III.2 and III.3. As a preliminary work, and in order to gain insight into the UJ-model in comparison with the stabilized jellium model, we study the stability of infinite uniform cylindrical UJ-wires.

We calculate the surface energy of the nanowires and the oscillations in the energy per unit length as it was made in our previous work describing Al, Na, and Cs nanowires within the stabilized jellium model Ogando et al. (2002); Puska et al. (2001). The results are shown in Fig. 1 as a function of the wire radius $`R`$. Here the radius is defined as the radius of the positive background charge in the SJ system with $`r_s`$ = 4.18 $`a_0`$ and the same amount of charge per unit length. In order to separate the energy oscillation from the average behavior, the so-called liquid-drop model Ziesche et al. (1993) is used. In this model the energy of the jellium system can be written as the sum of two terms – one proportional to the volume and the other proportional to the surface area. For the first term, the energy/volume ratio corresponds naturally to the homogeneous electron gasPerdew et al. (1991) with $`r_s`$ = 4.18 $`a_0`$. This view has been tested in clusters and nanowires Zabala et al. (1999); Puska et al. (2001); Ogando et al. (2002) and it describes correctly the mean energy, i.e. without the characteristic oscillations due to the quantum confinement. We fit the self-consistently calculated total energy per unit length to a liquid-drop model type function. Then, subtracting this smooth energy function from the total energy we get the pure energy oscillations, which are shown in the inset of Fig. 1. Note that there are radii for which the energy is at minimum. They correspond to wires which are more stable than wires with slightly different radii and higher energies. The first magic radii are at $`R=4.3`$, 7.3, 10.3, 13.6, 17.8, and 20.7 $`a_0`$. We use these radii for the initial uniform wires in the nanowire breaking simulation in Section III.3. The shell and super-shell structures studied in previous calculations Puska et al. (2001); Ogando et al. (2002) are also quite clear. In comparison with the energy oscillations of Na we observe that the beat positions are shifted to higher radii. The reason is that the UJ potential is softer at the surface than the stabilized jellium potential for Na Ogando et al. (2002).

### III.2  Periodic systems

Now we change the scheme and allow the wire to deform also in the axial direction. However, we impose periodic boundary conditions with the unit cell length $`L_{\mathrm{cell}}`$ along the wire axis. $`L_{\mathrm{cell}}`$ is thus the maximum perturbation wavelength in our calculation. From the liquid-drop model point of view, neglecting the small contribution of the curvature energy, the liquid wire attempts to achieve the shape which minimizes the surface, and thereby the total energy. Under this assumption an infinite periodic liquid wire is a uniform cylinder for lengths $`L_{\mathrm{cell}}<4.5R`$. For $`L_{\mathrm{cell}}>4.5R`$ it deforms trying to achieve the energetically most favorable state, an infinite chain of spheres. However, Kassubek et al. Kassubek et al. (2001) showed using a semi-classical model and perturbation theory that due to the discreteness of electronic structure the wires with magic radii remain uniform also at large $`L_{\mathrm{cell}}/R`$ ratios. With this result they argued that in the narrowing process of an infinite wire, when the radius is crossing an unstable zone before the next stable radius is achieved, the wire would spontaneously deform acquiring a wavy or deformed shape. We corroborate these results non-perturbatively using the UJ model as follows.

We choose a certain radius $`R`$ and solve for the UJ electronic (and positive charge density) structures imposing increasingly longer supercell lengths $`L_{\mathrm{cell}}`$ by increasing the number of electrons in the cell. Thereby we determine the critical supercell length (the wavelength of a perturbation) at which the wire starts deforming. For magic wires we find no wavy solutions, the wires remain uniform. For example, for $`R=7.3a_0`$ the wire is still uniform at $`L_{\mathrm{cell}}/R36`$. The wires corresponding to the radii at the maxima of the energy oscillations in Fig. 1 are the most unstable ones. These wires are uniform up to a critical value of $`L_{\mathrm{cell}}`$ but above it they spontaneously deform to a wavy or non-uniform density profile along the wire axis. As an example, Fig. 2 shows the behaviour of a wire with radius $`R=5.5`$ $`a_0`$ when the number of electrons in the unit cell is 6, 7, 8, and 10 and the unit cell lenght $`L_{\mathrm{cell}}`$ increases as 19.3, 22.5, 25.7 and 32.1 $`a_0`$, respectively. The unit cell with eight electrons correspond to a magic spherical cluster and that of ten electrons the pair of magic clusters of eight and two electrons. The critical values for the unstable radii of $`R=5.5`$, 8.6, 11.6 and 19 $`a_0`$ are $`L_{\mathrm{cell}}/R=`$4.1, 3.2, 4.2, and 4.8, respectively. I.e., we obtain values near the classical value of 4.5. At the unstable radius of $`R=15.5a_0`$ the wire is not deformed at least up to $`L_{\mathrm{cell}}/R=10`$ (the largest length we have calculated), probably due to the fact that this radius lies in a beat of the super-shell structure and it is actually relatively stable. We start all the calculations with a converged uniform potential profile along the wire axis (See Section III.1). In this way we do not “add any energy” to the system when initiating the calculation. Therefore, if the wire starts to deform in the iteration process the reason is the disappearance of the local energy minimum.

In addition, we narrow a stable uniform wire by increasing the length $`L_{\mathrm{cell}}`$ of the periodic cell and maintain the number of electrons constant. Each elongation step is solved self-consistently until convergence is reached. We observe that during the first steps the wire remains uniform but at some point, before breaking into isolated clusters, the wire spontaneously deforms. Thus, we confirm self-consistently and dynamically the hypothesis by Kassubek et al. Kassubek et al. (2001).

### III.3 Breaking of supported finite nanowires

In order to study the formation and evolution of nanoconstrictions between two supporting leads we follow the next procedure. First, we fix the number of electrons in the periodic supercell and solve self-consistently for the electronic structure of a uniform UJ wire having a stable magic radius. Then, the potential at both ends of the periodic cell is ”frozen”. This means that, although the Kohn-Sham equations are solved in the whole wire, in these regions the potential is not updated in the self-consistency process. The function of this ”frozen” part is to emulate the lead parts where ion rearrangement does not occur as efficiently as at the constriction. In our calculation, these leads serve as handles to grab the UJ and pull it. The rest of the wire, the UJ at the middle part of the supercell, is the place where the wire will stretch. A sketch of the configuration is shown in Fig. 3. A sharp change in the potential between the constriction and the leads turned out to cause difficulties in numerical calculations. Therefore, we smooth out the potential at the left edge using the form

$$F(zz_{\mathrm{edge}}^\mathrm{l})V_{\mathrm{frozen}}+F(z_{\mathrm{edge}}^\mathrm{l}z)V_{\mathrm{UJ}},$$
(8)
where F is a Fermi function with half-with of $`0.5a_0`$, $`V_{\mathrm{frozen}}`$ is the ”frozen” potential and $`V_{\mathrm{UJ}}`$ is the self-consistent UJ potential. For the right edge an analogous mixing is used. The main properties of the nanowire will not depend on the particular choice of this matching because the physical features are determined by the narrowest part of the constriction.

We perform simulations starting with radii between $`7.3a_0`$ and $`20.7a_0`$ and changing the number of electrons initially in the constriction. The elongation of the wire is made in steps of about one atomic unit, and always starting from the previous converged density, so that the grid spacing of the point mesh is increased to enlarge the cell. In order to overcome the interactions between the constriction and its periodic replica C.A. Stafford and Baeriswyl (2000), we choose the length of the lead part to be 6 Fermi wavelengths ($`\lambda _F=13.7a_0`$). Throughout the rest of the paper we will use $`\mathrm{\Delta }L`$ for the elongation; $`\mathrm{\Delta }L=0`$ for the first step.

In Fig. 4 we show snapshots of the electronic density for a wire with the starting radius of 10.7 $`a_0`$. The UJ-part corresponds to eight UJ-electrons placed initially in the neck region. Electrons are free to move inside or outside the leads, depending on the requirements of the self-consistent solution. However, there are always about eight electrons in the constriction. Although this is one of the smallest wires we have calculated, it shows all the main features observed when simulating also larger wires.

If the breaking of an UJ nanowire would happen as for fluid between the leads, the electron density should evolve forming a catenoid-shaped surface. Similar shapes (like hyperbolic Torres et al. (1994), parabolic Yannouleas et al. (1998), cosine Stafford et al. (1997), etc.) have been used before to model the nanoconstriction in simple free-electron or jellium simulations. The main results, when the comparison is possible, have been essentially the same irrespective of the actual shape. In Fig. 4(a) the electron density is shown after the elongation of $`\mathrm{\Delta }L=7.9a_0`$. The catenoid-like density profile appears as expected for a classical fluid. When we continue elongating the nanowire the shape of the electron density changes dramatically from the classical one. If the distance between the leads is short the electrons are strongly trapped at the narrowest part and they do not have much freedom in the rearrangement process. When the length of the constriction is large enough, the electrons have more space and freedom to achieve different types of energetically preferred shapes.

In Fig. 4(b) $`\mathrm{\Delta }L=19.8a_0`$ and the electrons in the constriction form a CDS. The electron density per unit length has two minima at both sides of the CDS and there are $`7.1`$ electrons between these narrowest cross sections. The embedded cluster reminds the closed-shell cluster of eight electrons, but there are some differences. There are not enough electrons and the symmetry is not exactly spherical. It seems that the $`p_z`$ orbital ($`z`$ along the cylinder axis) of the cluster has disappeared. We will analyse the structure in more detail below. Figure 4(c) shows the next consecutive elongation step with $`\mathrm{\Delta }L=20.8a_0`$. Note that the CDS disappears and a sudden change in the mean radius happens. In fact, the conductance changes simultaneously abruptly from $`3G_0`$ to $`1G_0`$ (See the inset in Fig. 5(a)). At this point it is also remarkable that the shape of the constriction is again far from the catenoid having a constant magic radius. Figure 4(d) is for $`\mathrm{\Delta }L=25.8a_0`$, the last step before the nanowire breaks. Again a CDS appears during the elongation from the third to the fourth snapshot. There are $`1.8`$ electrons between the two minimum cross sections at both sides of the CDS. This CDS can be interpreted as an embedded two-electron cluster. We observe that the radius of the constrictions is more or less constant with the same value as in the previous snapshot in Fig. 4(c).

At this point we want to focus on one characteristic property of UJ found when simulating the wire breaking: the UJ-matter deforms very easily. This ability to deform allows the formation of the cylinders of magic radii glued to the leads. The radius jumps from one magic radius to the next through an abrupt charge reorganization. The CDS’s of about two or eight electrons appear before the last charge reorganizations and the wire breaking. If there is enough UJ between the leads suspended long thin cylinders appear and in the last steps they alternate with chains of CDS’s producing a very extended elongation process. Here we want to underline that the CDS formation is a process different from the stability of a uniform cylindrical wire against the formation of a chain of spheres studied in the previous Section III.2. In that section, the quantum-mechanical shell structure may conserve the cylindrical structure which is not classically stable whereas now the quantum-mechanical shell-structure effect destroys the classical catenoid-type of solution producing a CDS in the constriction.

In Fig. 5 we show the conductance, the effective radius, and the elongation force as a function of the elongation for two different wires. The main figures correspond to an initial configuration with the radius of $`20.7a_0`$ and 60 electrons in the UJ-constriction. The insets display the results for a wire with an initial radius of 10.7 $`a_0`$ and eight electrons in the constriction. The electron density of the latter wire is plotted in Fig. 4 at certain elongation stages.

The conductance is calculated with the adiabatic and semiclassical approximation used by Brandbyge et al. Brandbyge et al. (1995). The constriction is divided into transversal slices. Then for each slice a uniform wire with the radial extent of the slice is built and the energy eigenvalues of the subband bottoms are calculated for this slice. The subband bottoms give effective potentials along the wire axis. If we look at the dependence of one of them on the position we see that it raises at the constriction due to the strong confinement (see Fig. 8). The electrons in this subband at the Fermi energy of the leads have to overcome this barrier in order to carry current. To evaluate the transmission probability of the electrons at the Fermi level through the barrier the semi-classical WKB formula is used.

The properties of the nanowires have been demonstrated to be dominated by the narrowest part of the constriction. Therefore we calculate an effective radius by evaluating the electron density per unit length at the middle of the wire. It is obtained with the value of the bulk electron density (corresponding to $`r_s=4.176a_0`$). Fig. 5(b) shows the effective radius as a function of the elongation of the wire. The plateaus or shoulders are in good coincidence with the infinite wire magic radii of $`10.3a_0`$, $`7.3a_0`$ and $`4.3a_0`$. For the larger wire shown in Fig. 5(b) also a small kink can be seen at $`\mathrm{\Delta }L=15.5a_0`$, which corresponds to the magic radius of $`13.6a_0`$. Wider magic radii do not appear because of the beat region of the supercell structure. In the inset at the end of the plateaus the effective radius increases when elongating the wire due to the CDS formation. The sudden decrease of the effective radius, accompanied by a step in the conductance, is due to the sharp charge rearrangements in the constriction.

The elongation force, shown in Fig 5(c) is evaluated as the negative derivative of the total energy with respect to the elongation. The rearrangement of the wire charge leads to discontinuous upward steps in the force, while if the radius changes smoothly the force draws a continuous buckling curve. At this point we want to point out the superiority of the UJ model in the force calculation over other jellium models Yannouleas and Landman (1997); Yannouleas et al. (1998); Zabala et al. (1999). In contrast with experiments Rubio et al. (1996); Rubio-Bollinger et al. (2001), the latter show a continuous behavior of the force without any steps. Moreover, for narrow constrictions positive values are obtained when the wire crosses an unstable zone. Note that in our model the force is always negative, as observed in the experiments Rubio et al. (1996); Rubio-Bollinger et al. (2001) and in atomistic simulations Landman (1998); Sørensen et al. (1998); Brandbyge et al. (1995); Nakamura et al. (1999). Fig. 5 shows clearly that the transport, geometrical and mechanical properties of the nanowires under elongation are related.

### III.4 Electronic cluster-derived structures

Let us now analyse more closely the CDS appearing in Fig. 4(b). In order to enlighten the origin of this structure we plot in Fig. 6(a) the single-particle energy spectrum of the wire. The extended zone scheme is used for clarity. The labels on the left of each branch represent the corresponding $`(|m|,n)`$ subbands for the infinite wire. In practice, $`n`$ is obtained by calculating the number of radial nodes at the cell boundaries (See Fig. 3). The branches have the characteristic parabolic shapes, but they show two different stages. In the lower part of the parabolic subbands the eigenvalues form flat plateaus without $`k_z`$ dispersion. These states correspond to wavefunctions localized at the leads and they vanish at the center of the constriction. Therefore the (0,2) and (2,1) subbands can not carry current through the constriction and they are closed channels. On the other hand, the states of the upper part of the (0,1) and (1,1) branches are extended along the whole wire and they form a continuous band (with the exception of small band gaps). The conductance of the wire is thus $`3G_0`$ due to the extended states of the $`(0,1)`$ and $`(1,1)`$ open channels at Fermi energy. This conclusion is in accordance with the value obtained with the WKB approach.

In Fig. 6(b) we plot the integrated local density of states (ILDOS) in the constriction for the band structure of Fig. 6(a). It is calculated by integrating the local density of states (LDOS) over the space between the two narrowest parts in the electron density in Fig. 4(b). The LDOS itself is obtained by substituting the discrete energy levels with Lorenzians of width (FWHM) of 0.4 eV and weighting them by the local probability amplitudes of the states in question. The ILDOS has two clear peaks, and when decomposing it we can see that the lower and the higher peak have the $`m=0`$ and $`m=1`$ character, respectively. The contribution of the $`m=2`$ states is negligible. The two ILDOS peaks can be fitted by two energy levels convoluted with the same Lorenzian as the eigenlevels in the LDOS calculation. The resulting resonance peaks are shown in Fig. 6(b) by dotted lines. The positions and the heights of these peaks have been fitted manually. The coincidence between the fit and the true ILDOS is remarkable. In the inset of Fig. 6(b) we plot the LDOS integrated between the leads for the electron density showed in Fig. 4(a) having no CDS. We observe that the ILDOS is much smoother and it is similar to the DOS of an infinite wire with delocalized states. The different $`m`$ contributions cannot be fitted by single resonance peaks as shown by the dotted peak for the $`m=0`$ and $`m=1`$ contributions. Moreover, the inset shows that the $`m`$-decomposed peaks are slightly asymmetric with a tail on the high energy side. These tails, which are not observable in the main figure in which the CDS appears, are due to the $`\sqrt{ϵ}`$-dependence of the subband peaks in the DOS for infinite wires.

The ILDOS analysis suggests that in the energy subbands or branches, at the transition points from states localized in the leads to states extended across the whole wire (see Fig. 6(a)), rather localized resonance states appear in the constriction. To clarify this point we plot selected states at the ILDOS peak energies in Fig. 7. Figs. 7(a) and 7(b) show clearly the localized character of the wavefunctions in the constriction at these energies. The state in Fig. 7(a) can be identified as the $`1s`$ orbital of an eight-electron cluster. The second well-localized state, Fig. 7(b), has $`m=1`$. Therefore it is doubly degenerate and it is identified as the $`p_{xy}`$ orbital. At about the energy of this $`p_{xy}`$ orbital a $`p_z`$ orbital (directed along the wire axis) should appear in the $`m=0`$ branch in order to complete the eight-electron cluster. However, we do not find such a state with a strong localization in the constriction. As shown in Fig. 7(c) the $`p_z`$ type states are much more delocalized than the $`p_{xy}`$ resonance states. The difference reflects the fact that due to the orientation the interaction of the cluster $`p_z`$ orbital with the lead states is much stronger than that of the $`p_{xy}`$ orbital. The absence of a well-localized $`p_z`$ orbital explains the clearly non-spherical shape of the embedded cluster in the electron density plot of Fig. 4(b), and also the finding that there are only 7.1 electrons in the constriction between the two narrowest cross sections.

Fig. 7(d) shows a well-localized state for the wire with eight UJ-electrons at the elongation of $`\mathrm{\Delta }L=25.8a_0`$ (Fig 4(d)). The state can be identified as the $`1s`$ orbital of a 2 electron cluster glued to the leads. There are 1.8 electrons between the two narrowest cross sections of the constriction supporting the assumption that this state is related to a two-electron cluster.

The states in Fig. 7 have always a wavy, non-decaying background. This is a characteristic of resonance states; truly localized states would decay exponentially. The wavy background corresponds to the wave function of the leads (plane wave) at the energy, which matches with that of the cluster state. To check this assumption we realize that the wavelength of the plane-wave background corresponds to the $`k_z`$ quantum number in the extended zone scheme: There are indeed two maxima in the modulus of the wave function per every Brillouin Zone unit of $`k_z`$ (see the labels of each wave function).

The existence of resonance states is related to the shape of the self-consistent potential having a small potential well in the nanoconstriction. To point out how this potential can admit a resonance state we show in Fig. 8 the effective potential for states with different $`(m,n)`$ quantum numbers, calculated within the adiabatic approximation for the wire with eight UJ-electrons and the elongation of $`\mathrm{\Delta }L=19.8a_0`$ (Fig. 4(b)). We see that electrons at the constriction feel the existence of a potential well. We plot the energies corresponding to the ILDOS peaks with dashed lines and note that they lie exactly in the potential wells, where the resonances situate. It is also evident that an occupied resonant $`p_z`$ state does not occur because its energy eigenvalue should be well above the effective potential of the $`(0,1)`$ branch and because the potential well of the $`(0,2)`$ branch is above the Fermi level. In addition, by the help of Fig. 8 we can explain the different parts of the electron energy bands in Fig. 6(a). The states with energies above the effective potential maxima are extended along the whole wire. These are the current-carrying states of each branch. The states below the potential minimum of the constriction are trapped in the leads, corresponding to flat plateaus in the lower part of the energy branches (Fig. 6(a)). Finally, between the potential maxima and the local minimum in the center we find resonant states which are enhanced at the constriction although they continue as plane waves in the leads.

## IV Conclusions

We studied the stability of nanowires and the nanowire breaking process performing self-consistent calculations within the ultimate jellium model. In the model, electrons and positive background charge acquire the optimal density minimizing the total energy. The model enables thus studies of shape-dependent properties of nanoscopic systems such as quantum dots or, as in the present work, quantum wires. The model advocates the idea that the electronic structure determines via the shell structure the geometry and ionic structure also in a partially confined system.

First we analysed the stability of infinite periodic quantum wires pointing out the ability of the electronic band structure to stabilize the nanowires at magic radii, i.e. any small deformation of the nanowire along the $`z`$ axis always increases the energy. At the unstable radii corresponding to maximum values of the energy oscillations the wire is uniform up to a critical value of the unit cell length. The critical values found are close to the classical value of $`L_{\mathrm{cell}}/R=4.5`$. Above this limit the local energy minimum disappears and a deformation of the wire lowers the total energy.

Then we investigated the elongation process of finite nanowires supported by leads. The elongation force, conductance and effective radius of the constriction were calculated simultaneously. The importance of the charge relaxation, in order to obtain results in agreement with experiments, was shown, e.g., in the case of the elongation force. The ability of the ultimate jellium (electron density) to acquire the optimal shape allows the formation of CDS’s showing the importance of electron states in the formation of these structures. The related resonance states and their origin was also shown. We found CDS’s that can be linked with the eight- and two-electron free-standing clusters.

###### Acknowledgements.

One of the authors (E.O.) acknowledges the Spanish Ministerio de Ciencia y Tecnología for financial support under the project PB98-0870-C02 and the Laboratory of Physics of the HUT for the kind hospitality. T.T. acknowledges financial support from the Vilho, Yrjö and Kalle Väisälä foundation. This work has also been supported by the Academy of Finland through its Centre of Excellence Program (2000-2005).


text:
# NUCLEOSYNTHESIS IN GAMMA-RAY BURST ACCRETION DISKS

## 1 INTRODUCTION

Growing evidence connects GRBs to the the birth of hyper-accreting black holes (Fryer, Woosley, & Hartmann, 1999), that is stellar mass black holes accreting matter from a disk at rates from $``$0.01 to 10 $`M_{\mathrm{}}\mathrm{sec}^1`$. Such models include the collapsar (Woosley, 1993; MacFadyen & Woosley, 1999), merging neutron stars and black holes (Eichler et al., 1989; Janka et al., 1999, 2001), supranovae (Vietri & Stella, 1998, 1999), merging helium cores and black holes (Zhang & Fryer, 2001), and merging white dwarfs and black holes (Fryer et al., 1999). For such high accretion rates the disk is optically thick, except to neutrinos, and very hot, consisting in its inner regions of a (viscous) mixture of neutrons and protons. At its inner boundary the disk connects to the black hole. Whatever processes accelerate the putative GRB-producing jet might therefore be expected to act upon some mixture of disk material and other background medium (e.g., the collapsing star in the collapsar model). Farther out, material will be lost from the disk in a vigorous wind (MacFadyen & Woosley, 1999; Narayan, Piran, & Kumar, 2002). In fact, Narayan et al. suggest that most of the disk will be lost to a wind except for those models and in those regions where neutrino losses dominate the energy budget.

It is thus of some consequence to know the composition of such disks. In the least case, the nucleosynthesis may be novel and could account for rare species in nature, such as the $`r`$-process. At most, the presence of free neutrons may affect the dynamics of the GRB jet (Derishev, Kocharovsky, & Kocharovsky, 1999; Fuller, Pruet, & Abazajian, 2000), the GRB neutrino signature (Bahcall & Mészáros, 2000), light curve (Pruet & Dalal, 2002), and afterglow (e.g. via a “pre-acceleration” mechanism similar to the one discussed by Beloborodov 2002). It is also of some consequence to know whether the disk wind consists of radioactive 56Ni, as is necessary if a visible supernova is to accompany the GRB. If the electron mole number, $`Y_e=\mathrm{\Sigma }(Z_iX_i/A_i)`$, is less than 0.485, the iron group will be dominated by 54,56Fe and other more neutron-rich species (Hartmann, Woosley, & El Eid, 1985) which will be incapable of illuminating the supernova. Since the jet itself is inefficient at heating sufficient matter to temperatures required for nuclear statistical equilibrium (T $`\text{ }>\text{ }5\times 10^9`$ K), supernovae seen in conjunction with GRBs (Galama et al., 1998; Bloom et al., 2002, for example) would be difficult to understand.

We have thus undertaken a survey of the nucleosynthesis that happens in the disks of rapidly accreting black holes. The work is greatly facilitated by the existence of numerical and semi-analytic solutions that yield the temperature-density structure and drift velocity (Popham, Woosley, & Fryer, 1999). These solutions have been verified in the case of the collapsar model by direct numerical simulation (MacFadyen & Woosley, 1999).

## 2 COMPUTATIONAL APPROACH

For a disk that is optically thin to neutrinos the neutron to proton ratio, $`n/p`$, is determined by the competition between electron and positron capture on nuclei and free nucleons. Since the baryon number per co-moving volume is conserved, it is convenient to work with the electron fraction $`Y_e`$ defined by

$$Y_e=\mathrm{\Sigma }Z_i(X_i/A_i)$$
(1)
where $`Z_i`$, $`X_i`$, and $`A_i`$ are the proton number, mass fraction, and atomic mass number (integer) of the species $`i`$. We shall refer to compositions with $`Y_e<`$ 0.5 as “neutron-rich”.

If we neglect lepton capture on bound nuclei and approximate the disk material as consisting of a mixture of free nucleons and $`\alpha `$-particles, the evolution of $`Y_e`$ with radius is given by

$$𝐮\left(Y_e\right)=V(r)\sqrt{\frac{12GM/rc^2}{1(V(r)/c)^2}}\frac{dY_e}{dr}=\lambda _{e^{}\mathrm{p}}\left(Y_e\frac{1X_n}{2}\right)+\lambda _{e^+\mathrm{n}}\left(1Y_e\frac{1X_n}{2}\right).$$
(2)
Here $`X_n`$ is the mass fraction of free nucleons, $`𝐮`$ is the 4-velocity of the flow, $`M`$ represents the black hole mass, and $`V(r)`$ is the radial drift velocity as measured in an inertial frame co-rotating with the disk. In the middle term in Eq. 2 we have assumed steady state conditions as well as cylindrical symmetry, and we have adopted the relatively simple Schwarzschild metric relevant for most of our calculations. Generally, relativistic effects do not make a big difference in our calculations of the electron fraction, though they are important for determing the structure of the disk.

In Eq. 2 $`\lambda _{e^{}\mathrm{p}}`$ and $`\lambda _{e^+\mathrm{n}}`$ are the rates for the processes

$`e^{}+\mathrm{p}`$ $``$ $`\mathrm{n}+\nu _e`$ (3)
$`e^++\mathrm{n}`$ $``$ $`\mathrm{p}+\overline{\nu }_e.`$ (4)
These rates are given by

$`\lambda _{e^{}\mathrm{p}}`$ $`=`$ $`K{\displaystyle _{\frac{\delta m}{m_e}}^{\mathrm{}}}w^2(w{\displaystyle \frac{\delta m}{m_e}})^2G_{}(1,w)S_{}(T,U_\mathrm{F})𝑑w`$ (5)
$`\lambda _{e^+\mathrm{n}}`$ $`=`$ $`K{\displaystyle _0^{\mathrm{}}}w^2(w+{\displaystyle \frac{\delta m}{m_e}})^2G_+(1,w)S_+(T,U_\mathrm{F})𝑑w.`$ (6)
Here $`K6.41410^4\mathrm{s}^1`$ determines the free neutron lifetime, $`\delta m1.293\mathrm{MeV}`$ is the neutron proton mass difference, and $`m_e`$ is the electron mass. The functions $`S_{}`$ and $`S_+`$ are the electron and positron distribution functions, while the $`G_\pm `$ are the Coulomb wave correction factors discussed in Fuller, Fowler, & Newman (1982).

The electron Fermi energy $`U_\mathrm{F}`$ is found by inverting the expression for the net electron number density,

$$n_e^{}n_e^+=\rho Y_eN_A=\frac{1}{\pi ^2}\left(\frac{m_ec}{\mathrm{}}\right)^3_0^{\mathrm{}}p^2(S_{}(T,U_\mathrm{F})S_+(T,U_\mathrm{F}))𝑑p$$
(7)
with $`p=\sqrt{w^21}`$.

We solve for the evolution of $`Y_e`$ by integrating eq. 2 for the results presented in Popham, Woosley, & Fryer (1999). Those authors assumed $`Y_e=1/2`$. Because we will show that in some cases $`Y_e1/2`$, and consequently that the electron Fermi energy and degeneracy pressure are substantially smaller than for $`Y_e=1/2`$, our results are not entirely self consistent, but should suffice given the very approximate nature of a one-dimensional calculation.

## 3 RESULTS

### 3.1 Scaling with Disk Viscosity and Accretion Rate

We present results for several different disk models. The parameters describing these models, as well as the electron fraction near the event horizon at $`r10^6\mathrm{cm}`$, are given in Table 1. In all cases we assume $`Y_e=1/2`$ at large radii. For simplicity we have adopted a constant black hole mass of 3 $`M_{\mathrm{}}`$ though the scaling relations for mass are obvious in Popham et al. All else being equal, a larger mass will give less electron capture.

Table 1 shows a clear trend with disk viscosity and accretion rate: electron capture in the inner disk becomes more pronounced with decreasing viscosity and increasing mass accretion rate. This arises because low viscosity flows inefficiently advect angular momentum outwards and are, for a given mass accretion rate, denser than high viscosity flows. Dense flows are electron degenerate, so that pair $`e^\pm `$ creation is suppressed and electron capture dominates over positron capture. Similarly a larger accretion rate also implies a denser disk - for a given value of $`\alpha `$ \- and more capture.

Before discussing in some detail the evolution of $`Y_e`$ in our calculations we note that the influence of black hole spin on the composition is also important. This is because the angular momemtum imparted to the black hole via the accreting matter will typically drive the Kerr parameter high. For example, a Kerr parameter of $`a0.9`$ is typical for collapsars. A larger Kerr parameter allows the disk to move to smaller radii before entering the event horizon, so that more electron capture and a smaller $`Y_e`$ will result. This is demonstrated by model G, which corresponds to a black hole with $`a=0.95`$. In this model $`n/p`$ near the event horizon is about ten times larger than $`n/p`$ for the same model except with $`a=0`$ (model B). The influence of the black hole spin on $`Y_e`$ is limited to regions quite near the event horizon. For $`r>\mathrm{\hspace{0.17em}10}^{6.5}\mathrm{cm}`$, $`Y_e`$ in model B is essentially the same as in model G. These considerations imply that as far as the composition is concerned, the main difference between flows around zero angular momentum holes and those around large angular momentum holes will be that the jet in the large angular momentum case will likely be more neutron rich. The composition in the bulk of the wind coming from the disk, however, will likely be similar in the two cases.

In Fig. 1 we show the evolution of $`Y_e`$ with radius for an accretion disk with $`\alpha =0.1`$ and a relatively low mass accretion rate (at least for standard collapsar models), $`\dot{M}=0.01M_{\mathrm{}}\mathrm{sec}^1`$. This flow is not dense enough to drive the electrons degenerate. Consequently $`Y_e`$ is governed by thermal $`e^\pm `$ capture. In this case the only asymmetry is the neutron-proton mass difference and $`Y_e`$ actually increases slightly owing to the threshold for the rate in Eq. 3, becoming greater than 0.5. The free nucleon mass fraction appearing in Fig. 1 is calculated from the simple estimate provided in Qian & Woosley (1996). Lepton capture on heavy nuclei is negligible for this disk.

Results for a hotter and denser flow, $`\dot{M}=0.1M_{\mathrm{}}\mathrm{sec}^1`$, $`\alpha =0.1`$, are shown in Fig. 5. The material at radii $`r10^{7.5}\mathrm{cm}`$ is mildly degenerate and the electron fraction is driven to $`Y_e0.44`$. The quantity $`\lambda _{e^{}\mathrm{p}}(r/V)`$ is the product of the electron capture rate and a rough measure of the dynamic timescale, or time left before the infalling material crosses the event horizon. Because $`Y_e`$ can not come to equilibrium unless $`\lambda _{e^{}\mathrm{p}}(r/V)1`$, $`Y_e`$ is not in equilibrium in this flow. Also, because $`\lambda _{e^{}\mathrm{p}}(r/V)1`$ when bound nuclei are present, weak processes on heavy nuclei are not important even in the extreme limit where the bound nucleons behave as free nucleons with respect to electron capture.

Figures 3 and 4 show the influence of viscosity on the composition of the flow. For $`\alpha =0.03`$ electron capture becomes important at $`r10^{7.6}\mathrm{cm}`$ and $`Y_e`$ is in close equilibrium with $`e^\pm `$ capture until $`r10^{6.5}\mathrm{cm}`$. There is at most $`1/10`$ of an electron capture per bound nucleus in the disk. For $`\alpha =0.01`$ the inflowing material becomes degenerate and neutron rich early on. Most of the free protons are locked into $`\alpha `$-particles until the $`\alpha `$-particles dissociate. Positron capture on the excess free neutrons results in a brief increase in $`Y_e`$, as seen in the bump at $`r10^{7.3}\mathrm{cm}`$ in Fig. 2. The final $`n/p`$ in this disk is very large, $`20`$.

The influence of the mass accretion rate on the composition is seen in Fig. 5, where we plot results for $`\dot{M}=1M_{\mathrm{}}\mathrm{sec}^1`$ and $`\alpha =0.1`$. This disk is quite similar to the $`\dot{M}=0.1M_{\mathrm{}}\mathrm{sec}^1`$, $`\alpha =0.3`$ case. Indeed, this is just what is expected from the scaling relations in Popham, Woosley, & Fryer (1999). For $`\dot{M}1M_{\mathrm{}}\mathrm{sec}^1`$ and larger, the assumption that the disk is optically thin to neutrinos begins to break down (MacFadyen & Woosley, 1999; Di Matteo , Perna, & Narayan, 2002). This influences both the structure of the disk, and, because neutrino capture becomes important, the evolution of $`Y_e`$. Because neutrino trapping likely influences the structure of the disk at the same level as the assumption that $`Y_e=1/2`$, we do not calculate $`Y_e`$ for disks calculated under the assumption of partial neutrino trapping (and $`Y_e=1/2`$).

### 3.2 Effects of Neutrino Capture

The effect of neutrino capture on the evolution of $`Y_e`$ can be addressed in an approximate way. In particular, we are interested in the rate $`\lambda _{\nu _e\mathrm{n}}`$ for the process

$$\nu _e+n\mathrm{p}+e^{}.$$
(8)
The rate for $`\overline{\nu }_e`$ capture is suppressed both because of the low proton number density and the absence of high energy positrons in the flow. The number of neutrinos captured per neutron is roughly $`\tau _{\nu _e,\mathrm{cap}}n_{\nu _e}/n_\mathrm{n}Y_e\tau _{\nu _e,\mathrm{cap}}n_{\nu _e}/n_\mathrm{p}`$, where $`\tau _{\nu _e,\mathrm{cap}}`$ is the optical depth for the process in Eq. 8, and $`n_{\nu _e}/n_\mathrm{n}`$ ($`n_{\nu _e}/n_\mathrm{p}`$) is the number of electron neutrinos produced per neutron (proton) in the flow. The ratio $`R`$ of the number of neutrino captures per neutron to the number of electron captures per proton is

$$R\frac{\lambda _{\nu _e}}{\lambda _{e^{}\mathrm{p}}}=\tau _{\nu _e,\mathrm{cap}}Y_e.$$
(9)
This equation implies that an optical depth of 1 to $`\nu _e`$ capture is roughly equivalent to a doubling of $`\lambda _{e^+\mathrm{n}}`$. We note here that while the the $`\nu _e`$’s produced in electron-capture reactions are not thermal, their average capture cross section is only approximately a factor of 2 higher than the average capture cross section calculated under the assumption of a thermal distribution for the electron neutrinos (at the electron Fermi energy and temperature) for conditions of interest in accretion disks.

It is difficult to quantitatively calculate the influence of neutrino capture on $`Y_e`$ in a post-processing step. This is both because of the complexities of neutrino transport, as well as because of the feedback between $`Y_e`$, the electron capture rate, and the disk dynamics. However, the following simple considerations argue that neutrino capture is unlikely to have a dramatic influence on $`Y_e`$.

To discuss some specific cases, consider first the $`\dot{M}=0.1M_{\mathrm{}}\mathrm{sec}^1`$, $`\alpha =0.1`$ disk. Di Matteo , Perna, & Narayan (2002) argue that neutrinos will have an absorptive optical depth of unity at $`r310^6\mathrm{cm}`$ for this disk. Because weak processes have essentially frozen out by this radius (see the $`\lambda _{e^{}\mathrm{p}}(r/V)`$ curve in Fig. 2), neutrino capture will have a small influence on $`Y_e`$.

Denser flows trap neutrinos more efficiently but are still expected to remain neutron rich. For example, for the $`\dot{M}=1M_{\mathrm{}}\mathrm{sec}^1`$, $`\alpha =0.1`$ disk, Di Matteo , Perna, & Narayan (2002) show that the neutrino absorptive optical depth will be greater than unity for $`r210^7\mathrm{cm}`$. As weak processes are rapid compared to the dynamic timescale in this disk, the result will be an increase in $`Y_e`$. To estimate the magnitude of this increase we adopt the simple procedure of artificially increasing the positron capture rate, $`\lambda _{e^+\mathrm{n}}5\lambda _{e^+\mathrm{n}}`$, which is likely an overestimate of the protonization rate. The increase in $`\lambda _{e^+\mathrm{n}}`$results in an electron fraction of $`Y_e0.166`$ near the event horizon. This is approximately 50% higher than the electron fraction calculated by neglecting the influence of $`\nu _e`$ capture, though still quite neutron rich.

To summarize, $`Y_e`$ in disks with high viscosity and modest accretion rates will not be affected by neutrino capture simply because neutrinos are not trapped or because weak processes are too slow. The electron fraction in disks with neutrino absorptive optical depths of a few will increase somewhat but will remain neutron rich owing to the sharp rise in $`\lambda _{e^{}\mathrm{p}}`$ with $`Y_e`$.

### 3.3 Transport to the Surface of the Disk

In order for a low electron fraction to have observable implications the electron fraction must not come to equilibrium at $`Y_e=1/2`$ as the material travels out of the plane of the disk. To estimate the evolution of the electron fraction in convective blobs moving out of the disk we parametrize the convective timescale by the turbulent convection speed (Narayan, Piran, & Kumar, 2002)

$$v_{\mathrm{turb}}\alpha v_K=6.3\times 10^8\alpha _1r_7^{1/2}\mathrm{cm}\mathrm{s}^1$$
(10)
where $`\alpha _1=\alpha /0.1`$ and $`r_7=r/10^7\mathrm{cm}`$. This implies a time to go one pressure scale height

$$\tau _{\mathrm{conv}}r/v_{\mathrm{turb}}16\alpha _1^1r_7^{1/2}\mathrm{ms}.$$
(11)
We also assume that the convective blob expands adiabatically. An estimate for the entropy per baryon in the disk is given by Qian & Woosley (1996):

$$s/k_b0.052\frac{T_{\mathrm{MeV}}^3}{\rho _{10}}+7.4+\mathrm{ln}\left(\frac{T_{\mathrm{MeV}}^{3/2}}{\rho _{10}}\right).$$
(12)
Here $`T_{\mathrm{MeV}}`$ is the temperature in MeV and $`\rho _{10}`$ is the density in units of $`10^{10}\mathrm{g}\mathrm{cm}^3`$. The first term on the right hand side of Eq. 12 is the contribution to the entropy from relativistic light particles ($`\gamma /e^\pm `$), while the next two represent the contribution from free nucleons. Eq. 12 is not appropriate when electrons are degenerate. For degenerate electrons the contribution of thermal $`e^\pm `$ pairs is small and the coefficient 0.052 should be closer to 0.02, representing just the entropy of the photons. However, changing the coefficient to 0.02 makes little difference for our purposes.

When relativistic particles dominate the entropy, the adiabat satisfies $`T^3\rho `$. When free nucleons dominate the entropy, the adiabat satisfies $`T^{3/2}\rho `$. The evolution of the electron fraction in a convective blob is crucially sensitive to the adiabat on which the blob travels.

To illustrate this, consider first a disk with $`\dot{M}=0.01`$, $`\alpha =0.1`$. In Fig. 6 we show the evolution of the electron fraction in matter originating from two different points in this disk and for different assumptions about $`\tau _{\mathrm{conv}}`$. The upper curves in Fig. 6 correspond to material originating from $`r=10^{6.5}\mathrm{cm}`$ in the disk. At this point the entropy is relatively high, $`s/k_b23`$ by the estimate in Eq. 12, and consequently the adiabat is closer to $`T^3\rho `$ than to $`T^{3/2}\rho `$. The temperature decreases quite slowly relative to the density, allowing for the formation of pairs which efficiently drive $`n/p`$ to equality. For large $`\tau _{\mathrm{conv}}`$ this process occurs more efficiently than for small $`\tau _{\mathrm{conv}}`$. The lower curves in Fig. 6 correspond to matter originating from $`r=10^7\mathrm{cm}`$ in the disk. At this point the entropy is $`17`$, so the adiabat is near $`T^{3/2}\rho `$. For low entropies, then, the temperature and density decrease at a comparable rate and pair formation is somewhat suppressed, keeping $`Y_e`$ low.

For flows where the electron fraction is high enough to drive the inner disk neutron rich, the entropy is generally dominated by non-relativistic particles, and the adiabat again preserves a low $`Y_e`$. This is illustrated in Fig. 7 where we plot the evolution of electron fraction in an adiabatic convective bubble for two different disks characterized by a low $`Y_e`$. In this figure we only show the evolution of a fluid element from a single initial radius ($`10^7\mathrm{cm}`$) because the entropy in these is disks never dominated by relativistic particles when $`Y_e`$ is low, so that the behavior shown in Fig. 7 is generic.

The above considerations outline the general features of how $`Y_e`$ evolves in material as it travels from the center to the surface of the disk. In flows where $`Y_e`$ is driven small, the degeneracy of the electrons implies that free nucleons dominate the entropy. A low $`Y_e`$ then, will remain low. By contrast, in flows where $`Y_e`$ is close to 1/2, the adiabat can be close enough to $`T^3\rho `$ for $`Y_e`$ to be driven to 1/2, and the neutron excess at the center of the disk is larger than the neutron excess in material finding its way out of the disk.

## 4 IMPLICATIONS

### 4.1 Radioactivity in the disk wind

The most obvious consequence of our calculations is that 56Ni will be absent from the winds of hyper-accreting black holes unless a) the accretion rate is low ($`\text{ }<\text{ }0.1M_{\mathrm{}}\mathrm{sec}^1`$), and b) the disk viscosity is high, $`\alpha \text{ }>\text{ }0.1`$. Interestingly, modern views regarding $`\alpha `$-disks and GRB models favor values close to these.

Typical accretion rates in the collapsar model are 0.05 - 0.1 $`M_{\mathrm{}}\mathrm{sec}^1`$ (MacFadyen & Woosley, 1999, their figures 5 and 10). The neutron excess will be smaller in Type II collapsars powered by fall back rather than direct black hole formation (MacFadyen, Woosley, & Heger, 2001). Considerable accretion into the hole and mass loss from the disk may continue, at a declining rate, even after the main GRB producing event ($``$20 s) is over (Zhang & Woosley, 2002). It thus seems likely that the collapsar model will be able to provide the 56Ni necessary to make the supernovae that accompany GRBs (though only if $`\alpha \text{ }>\text{ }0.1`$). This is also true of the slower accreting models like helium-star black hole mergers and black hole white dwarf mergers. However any wind from merging compact objects, or similar models like the supranova, will be neutron rich. Though perhaps of interest for nucleosynthesis, they will not produce 56Ni, at least during the black hole accretion epoch.

### 4.2 Neutron excess in GRB jets

For disks with high accretion rates, and certainly for merging neutron stars or black hole neutron star mergers, the matter near the event horizon will be very neutron rich. If this material pollutes the outgoing jet, the the GRB jet will itself, at least initially, contain free neutrons.

The dynamics of accelerating neutron-rich jets can differ dramatically from the dynamics of pure proton jets. Fuller, Pruet, & Abazajian (2000) showed that a high Lorentz factor fireball that is neutron rich can lead to two very distinct kinematic components, a slow neutron outflow and a fast proton outflow. This arises because the uncharged neutrons are weakly coupled to the radiation dominated plasma, and are accelerated principally via strong neutron-proton scatterings (Derishev, Kocharovsky, & Kocharovsky, 1999). Strong scatterings freeze out when they become slow compared to the dynamic timescale and at this point the neutrons coast. If this decoupling occurs while the jet is still accelerating, then the coulomb-coupled protons go on to have a larger Lorentz factor than the neutrons.

Roughly, dynamic neutron decoupling is only expected for fast jets. Here the precise meaning of “fast” depends, among other things, the fireball source size. For relativistic flows originating from compact objects, the dynamic timescale characterizing the acceleration of the flow is about 1 ms and the final Lorentz factor must be greater than $`300`$ in order for dynamic neutron decoupling to occur. For jets in the collapsar model, the timescale characterizing the acceleration is set by the surface of last interaction of the jet with the stellar envelope at $`10^{11}\mathrm{cm}`$. In this case the final Lorentz factor of the jet has to be $`3000`$ for dynamic neutron decoupling to occur. Because this is an impossibly high Lorentz factor for the collapsar model, dynamic neutron decoupling will not occur.

There may be a detectable neutrino signature of dynamic neutron decoupling. Neutron-proton collisions occurring during decoupling will generate pions. In turn these pions will decay and lead to the generation of neutrinos with energies of a few GeV in our reference frame. These neutrinos should be detectable at the rate of about one per year in next generation neutrino telescopes (Bahcall & Mészáros, 2000). In addition, there may also be a direct electromagnetic signature of neutron decoupling. For some bursts arising from external shocks, a slow decoupled neutron shell will decay and shock with the outer proton shell as the outer shell plows into the interstellar medium. This leads to a characteristic two-peaked structure in the burst (Pruet & Dalal, 2002).

Regardless of whether dynamic neutron decoupling occurs, neutrons in the jet may have observable implications for the GRB afterglow. This is particularly true when the shocking radius, i.e. the radius at which the shocks giving rise to the observed $`\gamma `$rays occur, is smaller than the length over which free neutrons decay. This condition is

$$r_{\mathrm{shock}}\gamma _2\tau _nc10^{15}\gamma _2\mathrm{cm},$$
(13)
where $`r_{\mathrm{shock}}`$ is the shocking radius, $`\gamma _2`$ is the Lorentz factor of the outflow in units of 100, and $`\tau _n1000\mathrm{sec}`$ is the free neutron lifetime. Eq. 13 is generally satisfied for bursts from internal shocks and bursts from external shocks in a very dense medium. When Eq. 13 holds, the neutrons can stream ahead of the slowing proton shell (the complement of the dynamic neutron decoupling discussed above) and deposit energy as they decay. Implications for the resulting afterglow in such a case have been discussed by Beloborodov (2002).

In the next section we discuss the r-process in winds and jets from accretion disks. However, we first note that GRB jets are characterized by interesting light element synthesis. Pruet, Guiles, & Fuller (2002) and Lemoine (2002) calculated the thermal synthesis of light elements in GRB-like outflows. They find that the final deuterium mass fraction can be of order 1%, some three orders of magnitude larger than the primordial yield. For kinematically well coupled flows, the nucleosynthesis depends sensitively on $`Y_e`$, with deuterium yield decreasing with increasing neutron excess. For flows in which dynamic neutron decoupling occurs, high energy neutron-$`\alpha `$ collisions will spall deuterons and can result in final deuteron mass fractions as high as 10% (Pruet, Guiles, & Fuller, 2002).

### 4.3 The r-process

There are two possible sites for the $`r`$-process here in those cases where $`Y_e`$ is low - in the jet and in the disk wind.

#### 4.3.1 In the jet

The jet has the merit of originating in the vicinity of the black hole where the neutron excess is likely to be greatest. If magnetic fields drive the heating and initial acceleration of the jet, then this neutron excess will likely be preserved. If neutrinos drive the outflow (MacFadyen & Woosley, 1999), $`Y_e`$ will be reset to some extent by the weak interactions. However, if the details of the outflow above the black hole are similar to spherically symmetric neutrino driven ultra-relativistic winds, then only $`\nu \overline{\nu }`$ annihilation will be important and the flow will remain neutron rich (Pruet, Fuller, & Cardall, 2001).

It is clear that the rapid expansion of the jet will be favorable to freezing out with a large abundance of free neutrons and this is conducive to the $`r`$-process (Hoffman, Woosley, & Qian, 1997). However, a pure nucleonic jet coming from the adiabatic expansion of a fireball with energy loading $`\eta =E_{\mathrm{internal}}/\rho c^2\text{ }>\text{ }100`$ is too much of a good thing when it comes to entropy and rapid expansion. The entropy was given in eq. 12 (Qian & Woosley, 1996)

$$s/k=\frac{11\pi ^2}{45}(\frac{kT}{\mathrm{}c})^3m_N/\rho $$
(14)
The internal energy (erg/gm) is

$$ϵ=\frac{11\pi ^2}{60}(\frac{kT}{\mathrm{}c})^41/\rho $$
(15)
and hence, for $`\eta =200`$ and $`kT`$ 2 MeV

$$s/k_b=\frac{4}{3}\frac{\eta c^2m_\mathrm{N}}{kT}10^5.$$
(16)
The time scale for the expansion, $`r/c`$, is less than a millisecond. The high entropy implies a low density at the temperature when $`n`$ and $`p`$ can start to recombine. Coupled with the rapid expansion, one might conclude that the jet will remain pure nucleons and light nuclei and that an $`r`$-process is impossible.

This conclusion would be erroneous on several counts. First, the jet may not always have such high energy loading. There may be cases where no gamma-ray burst is produced and the energy loading of the jet is much less. The entropy would then decrease and the expansion time scale would increase. Whether to call these “dirty fireballs” or simply an extension of the “disk wind” is simply a matter of taste (see below).

Of more relevance to GRB models, and particularly to those that involve massive stars, is that the jet will not escape the star without interaction. Along the walls of the jet, the Kelvin Helmholtz instability will occur (Zhang & Woosley, 2002). The jet will also encounter one or more shocks where it will be abruptly slowed. Most interesting is the “jet head” which moves through the star subrelativistically with a speed that varies but is of order c/3. At this head the jet mixes with stellar material and the mixture is swept backwards (in the moving frame) forming a cocoon. It typically takes 5 to 10 sec for the jet head to penetrate the roughly one solar radius of overlying Wolf-Rayet star in the collapsar model.

During this time, neutrons in the jet are mixed with heavy nuclei - He, O, Si, Ne, Mg and the like - which can serve as seeds for an $`r`$-process. The overall dynamics is likely to be quite complex and its study would requires a multi-dimensional relativistic simulation well beyond the scope of this paper. However, Zhang & Woosley (2002) find typical densities and time scales in the jet prior to break out are 10-4 to 10-1 g cm-3 and seconds. This is enough that many, if not all of the neutrons in the jet would capture.

Still, the overall yield in gamma-ray bursts is too small to account for the solar $`r`$-process. Assume there is one GRB for each 100 supernovae. Each burst has 1051 \- 1052 erg of relativistic ejecta (both jets) and a Lorentz factor that is at least 200 (Lithwick & Sari, 2001). The jet will be some mixture of entrained stellar material and jet, but assume that the initial jet is half or more of the relativistic ejecta. This gives an equivalent yield of about 10-7 $`M_{\mathrm{}}`$ per supernova. Probably the $`r`$-process will be at most 10% of this (the rest may be $`\alpha `$-particles and entrained nuclei). This then is several orders of magnitude less than required to make the solar r-process. However, GRB jets might still be important sources of the r-process in metal-deficient stars or of select rare nuclei in the sun.

#### 4.3.2 The $`r`$-process in the disk wind

A more significant $`r`$-process could come from the non-relativistic ejecta that comprise the disk wind, provided that the entropy is increased by magnetic processes to well above its value in the middle of the disk, as in for example, Daigne & Mochkovitch (2002). Assume that the accretion is 50% efficient, that is that half the material that enters the disk eventually enters the black hole and the other half is blown away (Narayan, Piran, & Kumar, 2002). A typical collapsar powered GRB involves the accretion of from one to several solar masses, so $``$1 $`M_{\mathrm{}}`$ is ejected from the disk. Most of this will not be neutron-rich enough or expand fast enough to make the $`r`$-process, but consider the implications if only 1% of the mass that is lost comes from the inner disk at times when the accretion rate is over 0.1 $`M_{\mathrm{}}\mathrm{sec}^1`$ and $`Y_e<0.4`$.

So close to the hole, significant entropy will be added by any acceleration process that produces a strong outflow. A crude estimate is $`s/k_bGMm_n/(rkT_{max})400r_6^1(2\mathrm{MeV}/T_{max})`$ (Qian & Woosley, 1996), where $`T_{max}`$ is the temperature where the energy deposition is the greatest. In fact the conditions in the inner disk for accretion rates $``$0.1 $`M_{\mathrm{}}\mathrm{sec}^1`$ are not so different - in terms of neutrino luminosity, neutron excess, temperature, and gravitational potential - from the neutrino-driven wind in an ordinary supernova, and one might expect a similar $`r`$-process (Woosley et al., 1994). The gravitational potential though, can in principle, be greater and a larger entropy and shorter time scale might be realized. Assume that the ejecta will be very rich in $`\alpha `$-particles and contain of order 10% $`r`$-process by mass. Putting the numbers together then gives an $`r`$-process synthesis equivalent to 10-5 $`M_{\mathrm{}}`$ per supernova, making this possibility well worth further investigation.

The authors acknowledge helpful correspondence with Weiqun Zhang regarding the conditions in collapsar jets. This research has been supported by NASA (NAG5-8128, NAG5-12036, and MIT-292701) and the DOE Program for Scientific Discovery through Advanced Computing (SciDAC; DE-FC02-01ER41176). A portion of this work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore Laboratory under contract W-7405-ENG-48.