text:
# 1 Introduction

## 1 Introduction

The most successful computational model for age-structured populations is by far the Penna model . It’s biological support comes from the mutation accumulation theory of senescence . In essence, this theory relates the evolution of senescence to the action of age-specific deleterious mutant alleles and the maintenance of some of these genes by the combined effects of mutation and selection pressures .

One of the reasons for the Penna model’s success relies on a particularly well-suited computational representation of a genome by means of a sequence of bits, the bit-string. When grouped into computer words these strings can be efficiently operated on by very fast logical and bitwise CPU operations. It has also proved to be flexible enough to be of value in a number of different problems in population dynamics, running from the catastrophic senescence of salmon to the origins of menopause (for a review we address the reader to references ; the first one contains also a fortran program to implement the model).

In order to avoid an exponential increase of the population, the Penna model makes use of the so-called Verhulst factor. It is a logistic-type term that introduces a mean-field random death probability, independent of the quality of the genome. Its usual expression is $`V(Pop)=Pop/Popmax`$, where $`Pop`$ is the total population at some time step and $`Popmax`$ is a parameter of the simulation, traditionally called the carrying capacity of the environment. Since really random deaths in nature can hardly play any significant role in population dynamics, this concept has already been criticized in the literature . Here we modify the model in such a way that individuals now die only for genetic reasons. Moreover, the population size is fixed, without fluctuations, as the biologists usually prefer.

## 2 Model with constant population

We start the simulation with $`N`$ individuals, each one represented by a chronological genome. Each of these genomes is reduced to the 32 bits of one computer word, where each bit represents a life-threatening inherited disease. The lifespan is divided into 32 intervals, each corresponding to one bit position. A zero bit means health; a bit set to one means that starting from that age interval until death, one additional disease is diminishing the individual’s health. $`T`$ such diseases, i.e. $`T`$ bits set to one in the bitstring from age zero to the current age, kill the individual. Those who survive up to age $`R`$, the minimum reproductive age, may get offspring, from then on. Each offspring differs in $`m`$ randomly selected bit positions from the parent. As usual, we consider only bad mutations which means that if a position is selected which has already its bit set to one, then it remains one in the offspring genome.

The reason why we underlined the word “may” in the above paragraph is because at that point we start to modify the standard model. Now as soon as an individual dies, either for reaching the age 32 or the limit $`T`$ for allowed diseases, a random individual is chosen, among all those with age $`R`$, to generate an offspring. In this way we keep the population size constant, without using the Verhulst factor already mentioned. This strategy can be interpreted as if, at each iteration, we have been chosing a finite sample of an infinite population to observe. In the next section we present our results, which are qualitatively the same as those obtained earlier by Sá Martins and Cebrat . However, in their paper they still keep a Verhulst factor acting over the newborns. For this reason their population size fluctuates a lot at each iteration.

## 3 Results

In the figures below we compare the results of the Penna model at constant population size with those obtained with the standard model. In order to have nearly the same population sizes, we have used for the standard model a carrying capacity $`Popmax=3,000,000`$ individuals, starting with $`P_0=50,000`$ individuals. With these parameters the population size of the standard model fluctuates around 360,000 individuals, after equilibration, while in the constant population model it is fixed to 300,000 since the beginning of the simulations. The rest of the parameters are the same for both: threshold for genetic deseases $`T=1`$, minimum reproduction age $`R=8`$, mutation rate of newborns $`m=1`$, total number of 1,000,000 timesteps with averages taken over the last 10,000.

In figure 1 we show the average fraction of indivuals per age. As obtained in , the overall concavity of the constant population curve is different from that of the standard one. Also, the maximum life span, the average age of an individual and the fraction of population with reproductive live (age above 8) are larger for the present model.

In figure 2 we present the fraction of defected genes (bits 1) for each locus in the genome (or age). As expected from the results of the previous figure, the fixation of deleterious mutations in the whole population genomes starts earlier with the standard model, causing a decrease of the maximum life span, when compared to the constant population model.

Figure 3 shows the mortality functions, $`q(a)`$, defined as:

$$q(a)=\mathrm{ln}\left[1\frac{D_a}{N_a}\right].$$

$`N_a`$ is the number of individuals with age $`a`$ and $`D_a`$ is the number of genetic deaths at age $`a`$ (which means, discounting the deaths provoked by the Verhulst factor in the standard model). It can be seen that the mortality above the minimum reproduction age increases faster in the standard model than in the present one. However, the standard model seems to give a better agreement with the Gompertz law , that predicts an exponential increase of the mortality above adult ages.

We have also measured the genetic diversity of the populations, computing the difference between each pair of genomes (Hamming distance). The diversity is higher for the constant population model, as also obtained in .

Real species, of course, can become extinct. Even within a constant population approach, by choosing an unrealistic set of parameters (for instance a high mutation rate), we observed mutational meltdown through a sudden explosion of the birth rate.

## 4 Conclusions

We removed from the Penna ageing model any random killing Verhulst factor and fixed the population size instead of allowing it to fluctuate around some constant value. We observe the same differences between this model and the standard one as those pointed out by Sá Martins and Cebrat: a longer life span, the fixation of deleterious alleles starting at older ages and a higher genetic diversity (not shown). However, we don’t have any fluctuation in the population size at each iteration, as they do, once we have a fixed number of individuals since the beginning of the simulations. The present model gives a mortality function that increases, for ages above the minimum reproduction age, in a slower way than that of the traditional Penna model. However, the standard model seems to give a better agreement with the Gompertz law of mortality.

Acknowledgements: To PMMH at ESPCI for the warm hospitality, to Sorin Tănase-Nicola for helping us with the computer facilities; to the Brazilian agencies FAPERJ and CNPq for financial support.


format: text

text:
# 1 Introduction

## 1 Introduction

We have shown earlier  that the spin of the “stable” mesons and baryons can be explained with the standing wave model . In  we have determined the rest mass of the $`\mu `$ mesons with the concepts used in the standing wave model and explained why m($`\mu ^\pm `$) is $``$ 3/4 $`\mathrm{m}(\pi ^\pm )`$. According to  the $`\mu ^\pm `$ mesons consist of a lattice of $`\nu _\mu `$ (respectively $`\overline{\nu }_\mu `$), $`\nu _e`$ and $`\overline{\nu }_e`$ neutrinos which remain from the cubic neutrino lattice of the $`\pi ^\pm `$ mesons after their decay, plus an electric charge. It has been argued that our explanation of the mass of the $`\mu `$ mesons cannot be correct because in our model the $`\mu `$ meson lattice has a diameter of 0.88$``$10<sup>-13</sup> cm, whereas the $`\mu `$ mesons are commonly said to be point particles. However, since in our model the $`\mu `$ mesons consist of neutrinos, but for the electric charge, and since neutrinos do not, in a first approximation, interact with electric charge or with mass, it will not be possible to establish the size of the $`\mu `$ meson lattice, i.e. of the $`\mu `$ mesons, through conventional scattering experiments. Therefore the $`\mu `$ mesons will *appear* to be point particles. We will now show that the neutrino lattice does not only determine the mass of the $`\mu ^\pm `$ mesons but also provides an explanation of the spin of the $`\mu `$ mesons which has, so far, not been explained.

## 2 The spin of the $`\mu ^\pm `$ mesons

The spin s = 1/2 or the intrinsic angular momentum $`\mathrm{}`$/2 of the $`\mu ^\pm `$ mesons can, theoretically, be the sum of the angular momentum vectors of all neutrino lattice oscillations of frequency $`\nu _i`$ in the $`\mu `$ mesons, plus the sum of all spin vectors of the $`n_i`$ neutrinos in the lattice, plus the spin vector of the electric charge which each $`\mu `$ meson carries. In a formula

$$j(\mu ^\pm )=\underset{i}{}j(\nu _i)+\underset{i}{}j(n_i)+j(e^\pm )0iN_\mu ,$$
(1)
where $`N_\mu `$ is the number of all neutrinos in the $`\mu `$ meson lattice, $`N_\mu =2.1410^9`$. This procedure is completely analogous to the way how the spin of the $`\pi ^\pm `$ mesons is determined, only that then a cubic neutrino lattice is considered consisting of $`\nu _\mu ,\overline{\nu }_\mu ,\nu _e`$ and $`\overline{\nu }_e`$ neutrinos with N = 2.854$`10^9`$ neutrinos.

The lattice oscillations in the $`\mu `$ mesons are longitudinal and hence do not have an angular momentum because for each oscillation $`\stackrel{}{r}\times \stackrel{}{p}=0`$. So the lattice oscillations do not contribute to the intrinsic angular momentum of the $`\mu `$ mesons, or $`_ij(\nu _i)=0`$. Each of the O($`10^9`$) neutrinos in the $`\mu `$ meson lattice has, however, an angular momentum $`\mathrm{}`$/2. The sum of the spin vectors of all neutrinos in the lattice of the $`\mu `$ mesons must be zero; otherwise the sum of the spin vectors of all neutrinos in the lattice plus the spin vector of the electric charge of a $`\mu ^\pm `$ meson could not be $`\mathrm{}`$/2, as it must be.

In order to show that the angular momentum vectors around a central axis of the lattice caused by the spin vectors of all neutrinos in the $`\mu `$ meson lattice is zero we have to consider this lattice in detail, Fig. 1. As we have learned in  the, say, $`\mu ^+`$ meson lattice is obtained from the cubic neutrino lattice of the $`\pi ^+`$ meson through the removal of all $`\nu _\mu `$ neutrinos. That means that in the $`\mu ^+`$ meson lattice are then N/4 = (2.854$`10^9`$)/4 vacancies at the location where originally $`\nu _\mu `$ neutrinos were. Each vacancy is surrounded in the xy plane by combinations of four electron neutrinos. In the z-direction there is on top as well as below each vacancy another electron neutrino. The same applies for the electron neutrinos around each antimuon neutrino, which are the centers of the cells of a $`\mu ^+`$ meson lattice. The cells of the $`\mu `$ meson lattice are octahedrons, two pyramids joint at their square base. As can be seen on Fig. 1 each lattice point in the upper right quadrant has at the opposite position in the lower left quadrant the *same* type of neutrino. The same applies in the upper left and lower right quadrant.

As is well-known it appears that only left-handed neutrinos and right-handed antineutrinos exist, at least as long as the neutrinos are massless. Suppose this also holds in the case of neutrinos with a rest mass. The angular momentum vectors originating from the spin of the neutrinos would then not cancel on a diagonal through the center of the lattice. The spin vectors on either side of the same length on the diagonal are then from the same type of neutrino (Fig. 1) and therefore point in the same direction and do not cancel. However, the polarization vector of the spin of a neutrino depends on the direction of the velocity of the neutrino, because the helicity H is given by

$$H=\stackrel{}{P}\stackrel{}{v}/Pv,$$
(2)
where $`\stackrel{}{P}`$ is the polarization vector. If only left-handed neutrinos (H = $``$$`\beta `$) and right-handed antineutrinos (H = $`+\beta `$) exist, the direction of $`\stackrel{}{P}`$ must be reversed if the direction of motion of the neutrinos during their oscillation in the $`\mu `$ meson lattice is, in the lower left quadrant of Fig. 1, opposite to the direction of motion in the upper right quadrant. This change of the direction of motion follows from the equation of motion for the displacements $`u_n`$ of the lattice points in Eq.(7) of

$$u_n=Ae^{i(\omega t+n\varphi )},$$
(3)
from which follows that $`\dot{u}_n=v_n=i\omega u_n`$. The frequencies are given by

$$\nu _n=\nu _0\varphi _n,$$
(4)
as in Eq.(19) of . Since n$`\varphi `$ = kx, $`\varphi _n`$ is proportional to x with x = n*a*, where *a* is the lattice constant. It follows that the direction of motion of the neutrinos in the upper right quadrant ($`\varphi >0`$) is opposite to the direction of motion of the neutrinos in the lower left ($`\varphi <0`$) quadrant. Consequently the angular momentum vectors around the center of the lattice caused by the spin of the neutrinos in the $`\mu `$ meson lattice are opposite and of the same magnitude, they cancel. The only point without an opposite is the point at the center of the lattice, at which there is no neutrino, so the center does not contribute a spin vector. The sum of the angular momentum vectors caused by the spin of all neutrinos in the $`\mu `$ meson lattice is zero, $`_ij(n_i)`$ = 0. Together with $`_ij(\nu _i)`$ = 0, as shown above, we arrive from Eq.(1) at

$$j(\mu ^\pm )=j(e^\pm ).$$
(5)
The spin s = 1/2 of the $`\mu ^\pm `$ mesons is caused exclusively by the spin of the electric charge that a $`\mu `$ meson carries. Crucial for this point is the absence of a neutrino at the center of the $`\mu `$ meson lattice.

The cubic lattice of the $`\pi ^\pm `$ mesons is easily recovered from the $`\mu `$ meson lattice by filling all vacancies with either $`\nu _\mu `$ (or $`\overline{\nu }_\mu `$) neutrinos. The $`\pi ^\pm `$ mesons do not have spin. As in the $`\mu `$ meson lattice all neutrino spin vectors around the center of the cubic $`\pi ^\pm `$ lattice cancel, but for the angular momentum $`\mathrm{}`$/2 of the center neutrino. This angular momentum is canceled by the spin of the electric charge which the $`\pi ^\pm `$ mesons carry, so s($`\pi ^\pm `$) = 0, as it must be. This explanation supercedes the explanation given in  for the spin of the particles of the neutrino branch which applies only for a static lattice.

## 3 Conclusions

If the $`\mu ^\pm `$ mesons consist of a lattice of $`\overline{\nu }_\mu `$ (respectively $`\nu _\mu `$), $`\nu _e`$ and $`\overline{\nu }_e`$ neutrinos, as is suggested by the decay of the $`\pi ^\pm `$ mesons, then it can be shown  that the theoretical mass of the $`\mu ^\pm `$ mesons is, within 1%, equal to the measured mass of the $`\mu ^\pm `$ mesons. Using the octahedronal lattice structure of the $`\mu ^\pm `$ mesons suggested by the determination of the mass of the $`\mu ^\pm `$ mesons it follows that the angular momentum vectors of all longitudinal lattice oscillations are zero and that the sum of the angular momentum vectors caused by the spin of the neutrinos of the lattice, taken around the center of the lattice, is also zero. The only contribution to the intrinsic angular momentum of the $`\mu `$ mesons comes from the spin of the electric charge which the $`\mu ^\pm `$ mesons carry. Consequently the spin of the $`\mu `$ mesons is s($`\mu ^\pm `$) = 1/2, as it must be. We note that the spin of the $`\mu `$ mesons can be explained, without any additional assumption, from the structure of the $`\mu `$ mesons which we have used for the explanation of the mass of the $`\mu `$ mesons.

Both the $`\pi ^\pm `$ mesons and the $`\mu ^\pm `$ mesons carry an electric charge. The $`\pi ^\pm `$ mesons do not have spin, whereas the $`\mu ^\pm `$ mesons have spin 1/2. The presence or absence of a neutrino at the center of the lattice makes the difference. The spin of the electric charge in $`\pi ^\pm `$ is canceled by the spin of the central neutrino, whereas the spin of the electric charge in $`\mu ^\pm `$ remains because there is no central neutrino to cancel the spin of the electric charge.

REFERENCES

 E.L.Koschmieder, arXiv: physics/0301060 (2003), Hadr.J. (to appear).

 E.L.Koschmieder, arXiv: physics/0211100 (2002), Chaos, Solitons and
Fractals 18,1129 (2003).

 E.L.Koschmieder, arXiv: physics/0110005 (2001).


text:
# Have the Elusive Progenitors of Supernovae Type Ia Been Discovered?

## 1 Introduction

The recent detection of H$`\alpha `$ emission in the spectrum of the supernova Type Ia (SN Ia) SN 2002ic (Hamuy et al. 2003) is a landmark discovery. While there is very little doubt that SNe Ia represent the thermonuclear disruption of mass accreting white dwarfs (WDs), the precise nature of the progenitor systems remains uncertain (Branch et al. 1995; Livio 2001). Given that SNe Ia are the tool of choice for confirming the acceleration of cosmic expansion (Riess et al. 1998; Perlmutter et al. 1999), the importance of identifying the progenitors cannot be overemphasized. The two main scenarios that have been proposed involve either the merger of two white dwarfs (the double-degenerate scenario; Iben & Tutukov 1984; Webbink 1984), or a single white dwarf accreting from a normal companion (the single-degenerate scenario; Whelan & Iben 1973; Nomoto 1982). Recently it has been argued theoretically, that single-degenerate progenitors are favored (even though it is very difficult for hydrogen-accreting WDs to reach the Chandrasekhar limit; Piersanti et al. 2000), and that double WD mergers may lead to accretion-induced collapses rather than to SNe Ia (Livio 2001; Nomoto et al. 2000). The tentative discovery (if confirmed) of an enhanced SN Ia rate near jets in active galactic nuclei (Livio, Riess, & Sparks 2002; Capetti 2002) appears to support this conclusion. Nevertheless, until SN 2002ic the “smoking gun”—the presence of hydrogen in the spectrum—was missing. The clear detection of a broad (FWHM $`1800`$ km s<sup>-1</sup>) H$`\alpha `$ component in SN 2002ic appears on the face of it to demonstrate that at least some SNe Ia result from single-degenerate progenitors. In the present letter we show that this conclusion may be *premature*.

## 2 Why Now?

One of the key questions posed by the observations of Hamuy et al. (2003) is: Why was hydrogen not detected before? This becomes particularly puzzling when we realize that there exist about 100 spectra of SNe Ia in which a signature of the strength of that seen in SN 2002ic would have been detected (T. Matheson, private communication), had it been there. In fact, Hamuy et al. noted that the amount of shock-heated circumstellar material needed to produce the observations of SN 2002ic is totally unexpected for a SN Ia. Accordingly, they suggested that the progenitor system was a binary consisting of a C/O white dwarf and a *massive* (3–7 M) asymptotic giant branch (AGB) star. The presence of the latter was necessitated by the need to have an integrated circumstellar mass of at least a few solar masses.

The main problem with this scenario is that one would expect to observe *a range* of strengths of H$`\alpha `$ lines in SNe Ia, depending on the amount of circumstellar material (in turn, determined primarily by the mass of the AGB star), rather than *detecting a relatively strong line in one case only* (it is also hard to believe that this is the first progenitor system containing an AGB star).

We propose instead that the total absence of H$`\alpha `$ lines in all the pre-SN 2002ic SNe Ia observed to date argues that SN 2002ic represents rather rare circumstances, and *not* a white dwarf accreting from the wind of an AGB star.

## 3 A Supernova Ia in a Common Envelope?

All the evolutionary scenarios leading to the formation of close double white dwarf systems involve a stage in which an AGB star fills its Roche lobe and transfers mass onto a white dwarf companion (e.g. Yungelson & Livio 2000). Under these conditions, the mass transfer process is unstable, and the system evolves rapidly into a common envelope (CE) configuration, inside which the white dwarf and the core of the AGB star spiral-in (e.g. Rasio & Livio 1996; Taam & Sandquist 2000). Typically, the CE phase lasts a few hundred to a few thousand years, and results in the ejection of the envelope and the emergence of a double white dwarf system (e.g. Sandquist et al. 1998; Taam & Sandquist 2000 and references therein). I propose that SN 2002ic represents one of those rare cases in which the explosion occurs *during* (or immediately following) the CE phase, and in which some part of the envelope has not been previously ejected. This raises two immediate questions: (i) Is this possible at all? and (ii) Does this support a single-degenerate or a double-degenerate scenario?

For the white dwarf to actually reach the Chandrasekhar mass via accretion of hydrogen-rich material during the CE phase is extraordinarily unlikely. Steady burning occurs for a narrow range of accretion rates of order (Paczyński & Żytkow 1978; Nomoto, Nariai & Sugimoto 1979; the limits are determined: at the low end by the requirement that the pressure at the time of ignition be sufficiently low to prevent a shell flash, and at the high end by the accretor expanding to supergiant dimensions)

$$0.4\dot{M}_{\mathrm{RG}}\dot{M}\dot{M}_{\mathrm{RG}}.$$
(1)
Here $`\dot{M}_{\mathrm{RG}}`$ is the rate at which the white dwarf expands to giant dimensions and is given by

$$\dot{M}_{\mathrm{RG}}8.5\times 10^7(M_{\mathrm{WD}}/\mathrm{M}_{}0.52)\mathrm{M}_{}/\mathrm{yr}.$$
(2)
Even assuming that the accretion rate could be regulated to the rate given by equation (1) \[most likely it would settle on the Eddington rate of $`\dot{M}_{\mathrm{EDD}}1.7\times 10^5(R_{\mathrm{WD}}/10^9\mathrm{cm})\mathrm{M}_{}/\mathrm{yr}`$ at which mass would not be retained\], the WD would increase in mass by at most $``$0.001 M during the CE phase. This would require the WD to be within 0.001 M of the Chandrasekhar mass upon entering the CE—a very unlikely situation, even taking into account the rarity of H$`\alpha `$ detection (e.g. only 2 out of a sample of 130 WDs were found to have masses higher than 1.2 M; Bergeron, Saffer, & Liebert 1992; although see Hachisu & Kato 1999).

A second possibility is that the WD spirals-in all the way to the center and merges with the AGB star’s core. Interestingly, a scenario for SNe of similar type was suggested almost 30 years ago by Sparks & Stecher (1974), but has long since been discarded due to the absence of hydrogen in the spectra. What I propose here is that the spiraling-in process unbinds most, but not all of the envelope, so that coalescence becomes inevitable. At the time of merger, most of the envelope will be at a distance of

$$d3\times 10^{15}\left(\frac{V}{10\mathrm{km}\mathrm{s}^1}\right)\left(\frac{\tau _{\mathrm{CE}}}{100\mathrm{yr}}\right)\mathrm{cm}.$$
(3)
from the core. Here $`V`$ is the ejection velocity and $`\tau _{\mathrm{CE}}`$ is the duration of the CE phase. The condition for a merger to occur (as opposed to ejection of the entire envelope and the formation of a binary WD system) is given by the requirement that the binding energy of the CE be larger than the gravitational energy available from orbital shrinkage (Livio 1996; deKool 1990)

$$\frac{M_{\mathrm{AGB}}(M_{\mathrm{AGB}}M_C)}{\lambda a_0r_L}>\alpha _{\mathrm{CE}}\left(\frac{M_CM_{\mathrm{WD}}}{2R_C}\frac{M_{\mathrm{AGB}}M_{\mathrm{WD}}}{2a_0}\right).$$
(4)
Here $`a_0`$ is the initial separation, $`r_L`$ is the Roche lobe radius of the AGB star (in units of the separation), $`M_C`$ and $`R_C`$ are the mass and radius of the core, respectively, $`\alpha _{\mathrm{CE}}`$ is the CE efficiency parameter (Livio & Soker 1988; Iben & Tutukov 1984), and $`\lambda 0.5`$ depends on the stellar density profile. The value of $`\alpha _{\mathrm{CE}}`$ is not known even to within a factor 10 (e.g. Livio 1996). However, for reasonable values ($`\alpha _{\mathrm{CE}}0.1`$–1) condition (4) requires relatively massive AGB stars \[since the condition can be approximated as $`(M_{\mathrm{AGB}}/M_{\mathrm{WD}})^21/8\alpha _{\mathrm{CE}}(a_0/R_C)`$; and $`a_0/R_C10^4`$\] and can be expected to be satisfied only in a fraction of a percent of all systems (e.g. Yungelson & Livio 1998). The observed H$`\alpha `$ emission would result from the interaction of the explosion with the previously-ejected envelope. This would be consistent with the rarity of H$`\alpha `$ detections. Most importantly, however, if this scenario is correct, the H$`\alpha `$ detection by Hamuy et al. *results from a double-degenerate scenario*!

## 4 Conclusions

One might have thought that the detection of hydrogen in the spectrum of a SN Ia would have finally revealed the elusive progenitor to be a single-degenerate system. In this Letter we suggest that this may not be the case. Paradoxically, the H$`\alpha `$ detection could result from a double-degenerate scenario! To be sure, the actual result of the merger process remains as uncertain as ever, and it may lead to an accretion-induced collapse rather than to a SN Ia. Other exotic possibilities, such as the explosion of the core of an AGB star (“type 1.5” event; Iben & Renzini 1983) may exist (as already suggested by Hamuy et al. 2003). However, the latter would require some other mechanism to place (at least a part of) the envelope at $`10^{15}`$ cm. Future, more sensitive, observations will reveal whether the detection of H$`\alpha `$ is a very rare, but relatively clear event, or whether a range of line strengths is detected. The latter case would clearly support a single-degenerate interpretation.

We would like to thank David Branch and Tom Matheson for helpful discussions.


text:
# Untitled Document

Erratum: ”Oscillations in Arcturus from WIRE photometry” (ApJ, 591, L151 )

The value of the mode lifetime given in the first paragraph of the Discussion should be divided by 2$`\pi `$ and is therefore 2.0 days rather than 13 days.


text:
# Spin-dependent Bohm trajectories for hydrogen eigenstates

## 1 Introduction

In David Bohm’s original causal interpretation of quantum mechanics , the motion of a quantum mechanical particle is determined by its wavefunction $`\psi `$, which acts as a guidance wave . If the wavefunction is written as

$$\psi (𝐱,t)=R(𝐱,t)e^{iS(𝐱,t)/\mathrm{}},$$
(1)
where and $`R`$ and $`S`$ are real-valued, then the trajectory of the particle is determined by the relation

$$𝐩=S.$$
(2)
The Schrödinger current associated with $`\psi `$ is given by

$$𝐣=\frac{1}{m}\rho 𝐩$$
(3)
where $`\rho =\psi ^{}\psi =R^2`$. Comprehensive discussions of this causal interpretation of quantum mechanics can be found in  and .

However, as Holland  has pointed out, Eqs. (2) and (3) are relevant only to spin-0 particles. For particles with spin, these equations are inconsistent if the theory is to be ultimately embedded in a relativistic theory. The condition of Lorentz covariance on the law of motion implies that the momentum of a particle with spin s must be given by

$$𝐩=S+\mathrm{log}\rho \times 𝐬.$$
(4)
In a number of earlier works, including , the current vector associated with Eq. (4),

$$𝐣=\frac{1}{m}\rho S+\frac{1}{m}\rho \times 𝐬,$$
(5)
was referred to as the Pauli current, the nonrelativistic limit of the Dirac current, as opposed to Eq. (3), the nonrelativistic limit of the Gordon current. (In these papers, it was claimed that consistency with Dirac theory requires that Schrödinger theory be regarded as describing an electron in an eigenstate of spin.) The spin-dependent term was also discussed in  but only in the context of the Pauli equation and not the Schrödinger equation.

The momentum defined in Eq. (2) predicts that electrons in an eigenstate are stationary since $`S=0`$. This is a counterintuitive result of Bohm’s original theory, but one which no longer persists when the extra term in (4) is taken into account. For example, consider an electron in the the $`1s`$ ground eigenstate of hydrogen,

$$\psi _{100}=\frac{1}{\sqrt{\pi a^3}}e^{r/a},$$
(6)
where $`a=\mathrm{}^2/(me^2)`$ is the Bohr radius. Also assume that the electron is in a definite spin eigenstate: Without loss of generality, let its spin vector be given by $`𝐬=\frac{\mathrm{}}{2}𝐤`$. (We shall justify this assumption below.) Holland  showed that the extra term $`\mathrm{log}\rho \times 𝐬`$ implies that the polar coordinates $`r`$ and $`\theta `$ are constant, with angle $`\varphi `$ evolving in time as follows:

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{mar}.$$
(7)
All points on a sphere of radius $`r`$ orbit the $`z`$-axis at the same rate. If $`r=a`$, the angular frequency is on the order of $`10^{16}\text{s}^1`$, cf. Eq. (32).

In Section 2, we examine the effects of the extra term in (4) for an electron in several eigenstates of the hydrogen atom. Because the momentum equation (4) now involves spin, a complete description of the electron in the atom – provided by an appropriate wavefunction – will have to involve both spatial as well as spin information. Let us denote the complete wavefunction of the electron by $`\mathrm{\Psi }(𝐱,𝐬,t)`$, where $`𝐬`$ denotes appropriate spin coordinates. Since the hamiltonian describing the evolution of $`\mathrm{\Psi }`$ is the simple spin-independent hydrogen atom hamiltonian $`\widehat{H}_0`$, we may write $`\mathrm{\Psi }`$ as the tensor product $`\psi (𝐱,t)\zeta (𝐬)`$, where $`\psi (𝐱,t)`$ is a solution of the time-dependent Schrödinger equation,

$$i\mathrm{}\frac{\psi }{t}=\widehat{H}_0\psi ,$$
(8)
and $`\zeta (𝐬)`$ is an eigenfunction of the commuting spin operators $`\widehat{S}^2`$ and $`\widehat{S}_z`$, with $`\widehat{S}^2\zeta =\frac{3}{4}\mathrm{}\zeta `$ and $`\widehat{S}_z\zeta =\frac{1}{2}\mathrm{}\zeta `$. Thus $`\zeta `$ defines the “alpha” or “spin up” state to which corresponds the spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$. As such, the remainder of our discussion can focus on the evolution of the spatial portion of the wavefunction $`\psi (𝐱,t)`$. In the parlance of “state preparation,” one can view this construction as preparing the electon with constant spin vector $`𝐬`$ and initial wavefunction $`\psi (𝐱,0)`$. Given an initial position $`𝐱_0`$ of the electron, Eq. (4) will then determine its initial momentum $`𝐩_0`$, from which its causal trajectory then evolves.

In Section 3, we examine the trajectories of an electron, as dictated by Eq. (4), where $`\psi `$ is a linear combination of $`1s`$ and $`2p_0`$ eigenstates evolving in time under the hydrogen atom hamiltonian, cf. Eq. (8). Once again because of the spin-independence of $`\widehat{H}_0`$, we may assume that the spin vector of the elctron is constant, i.e., $`𝐬=\frac{\mathrm{}}{2}𝐤`$, and thereby focus on the time evolution of the spatial wavefunction $`\psi (𝐱,t)`$. This simple model was chosen to determine the major qualitative features of Bohmian trajectories for a $`1s`$-$`2p`$ transition induced by an oscillating electric field. We have analyzed trajectories corresponding to the time-dependent wavefunction associated with the transition hamiltonian and shall report the results elsewhere.

## 2 Spin-dependent trajectories of electrons in hydrogen eigenstates

### 2.1 Qualitative features

Following the discussion at the end of the previous section, we consider an electron with spin vector $`𝐬`$ that begins in a hydrogenic eigenstate, i.e., $`\psi (𝐱,0)=\psi _{nlm}(𝐱)`$. The time evolution of the spatial wavefunction is simply

$$\psi (𝐱,t)=\psi _{nlm}(𝐱)e^{iE_nt/\mathrm{}}.$$
(9)
Comparing Eqs. (9) and (1), we see that for a real eigenstate, $`S=0`$ so that the momentum in Eq. (4) is given by

$$𝐩=\mathrm{log}\rho \times 𝐬.$$
(10)
Some simple qualitative information about the electron trajectories is readily found from this equation. First, the vector $`\mathrm{log}\rho `$ points in the direction of the steepest increase in $`\mathrm{log}\rho `$, hence in $`\rho =|\psi |^2`$. Because of the cross product, the momentum vector p is perpendicular to this direction. In other words, the trajectories of the electron lie on level surfaces of $`|\psi |^2`$. However, p is also perpendicular to the direction of the spin, assumed to lie along the $`z`$-axis in this discussion. This implies that z is constant for these Bohm trajectories.

From the above analysis, the shape of the Bohm trajectories may be found by computing level surfaces of $`|\psi |^2`$ – or simply $`\psi `$ for real-valued eigenfunctions – and then finding the intersections of these surfaces with planes of constant $`z`$. An electron in the $`1s`$ state of hydrogen, cf. Eq. (6), must therefore execute a circular orbit about the $`z`$-axis. However, the angular velocity of this orbit, cf. Eq. (7), cannot be determined from this analysis.

For the $`2s`$ case, with quantum numbers $`(n,l,m)=(2,0,0)`$, the wavefunction is given by

$$\psi _{200}=\frac{1}{\sqrt{32\pi a^5}}(1\frac{r}{2a})e^{r/2a}.$$
(11)
The level surfaces of $`\psi _{200}`$ are spheres whose intersection with planes of constant $`z`$ are circles (with constant $`\theta `$ values). The electron again travels about the $`z`$-axis in a circular orbit.

In the $`2p_0`$ case, $`(n,l,m)=(2,1,0)`$, the wavefunction is given by

$$\psi _{210}=\frac{1}{\sqrt{32\pi a^5}}re^{r/2a}\mathrm{cos}\theta .$$
(12)
The condition that both $`\psi `$ and $`z`$ be constant is satisfied only if both $`r`$ and $`\theta `$ are constant, once again yielding circular orbits about the $`z`$-axis.

In the other $`2p`$ cases, $`(n,l,m)=(2,1,\pm 1)`$, the wavefunctions are given by

$$\psi _{21(\pm 1)}=\frac{1}{\sqrt{32\pi a^5}}re^{r/2a}\mathrm{sin}\theta e^{\pm i\varphi }.$$
(13)
The condition $`\rho =|\psi |^2=\text{constant}`$ yields a relation $`r=r(\theta )`$ defined implicitly by

$$\mathrm{sin}\theta =\frac{K}{r}e^{r/2a},$$
(14)
where $`K`$ is a constant. Since $`z=r\mathrm{cos}\theta `$ is also constant, it follows that both $`r`$ and $`\theta `$ are constants of motion so that there is circular motion about the $`z`$-axis.

### 2.2 More detailed dynamical descriptions of the trajectories

For more complicated hydrogen eigenstates (see below), it may not be as simple to find closed form expressions for level sets of $`|\psi |^2`$ so that the method of qualitative analysis outlined above may be difficult if not impossible. As well, it is desirable to extract quantitative information such as the angular velocity $`d\varphi /dt`$ of the circular orbits deduced earlier. We therefore analyze the differential equations of motion defined by Eq. (10).

The gradient term from Eq. (10) is given by

$$\mathrm{log}\rho =\mathrm{log}\psi ^{}\psi =2\text{Re}\left(\frac{(\psi )\psi ^{}}{\psi \psi ^{}}\right).$$
(15)
For real wavefunctions, this simplifies to

$$\mathrm{log}\rho =2\frac{\psi }{\psi }.$$
(16)
In spherical polar coordinates $`(r,\theta ,\varphi )`$, the spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$ is given by $`𝐬=\frac{\mathrm{}}{2}(\mathrm{cos}\theta \widehat{r}\mathrm{sin}\theta \widehat{\theta })`$. It is convenient to compute the cross product in the (right-handed) spherical polar coordinate system:

$$𝐀\times 𝐁=\left|\begin{array}{ccc}\widehat{\theta }& \widehat{\varphi }& \widehat{r}\\ A_\theta & A_\varphi & A_r\\ B_\theta & B_\varphi & B_r\end{array}\right|.$$
(17)
For the $`1s`$ ground state defined in Eq. (6), $`\psi =\frac{1}{a}\psi \widehat{r}`$. Thus

$$𝐩=\frac{\mathrm{}}{a}\mathrm{sin}\theta \widehat{\varphi }.$$
(18)
Since $`p_r=p_\theta =0`$, it follows that $`r`$ and $`\theta `$ are constant, implying that $`z`$ is constant, i.e., circular orbits about the $`z`$-axis. Holland’s result in Eq. (7) follows.

For the $`2s`$ wavefunction defined in Eq. (11),

$$\frac{\psi _{200}}{\psi _{200}}=\frac{1}{2a}\left[\frac{1}{1\frac{r}{2a}}+1\right]\widehat{r},$$
(19)
so that

$$𝐩=\frac{\mathrm{}}{2a}\mathrm{sin}\theta \left[\frac{1}{1\frac{r}{2a}}+1\right]\widehat{\varphi }.$$
(20)
Once again $`r`$ and $`\theta `$ are constant. From the relation $`d\varphi /dt=p_\varphi /mr\mathrm{sin}\theta `$, we have

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{2mar}\left[\frac{1}{\frac{r}{2a}1}1\right].$$
(21)
Note that the pole at $`r=2a`$ coincides with the zero of the $`2s`$ wavefunction, implying that the probability of finding the electron at $`r=2a`$ is zero. Also note that (i) $`\dot{\varphi }>0`$ for $`0<r<2a`$, (ii) $`\dot{\varphi }<0`$ for $`2a<r<3a`$, (iii) $`\dot{\varphi }=0`$ for $`r=3a`$ and (iv) $`\dot{\varphi }>0`$ for $`r>3a`$. For $`r=a`$, the angular velocity $`\dot{\varphi }`$ is equal to that of the $`1s`$ ground state, cf. Eq. (7).

For the $`2p_0`$ state defined in Eq. (12),

$$\frac{\psi _{210}}{\psi _{210}}=\frac{1}{r}(1\frac{r}{2a})\widehat{r}\frac{\mathrm{sin}\theta }{r\mathrm{cos}\theta }\widehat{\theta }.$$
(22)
From Eq. (10),

$$𝐩=\frac{\mathrm{}\mathrm{sin}\theta }{2a}\widehat{\varphi },$$
(23)
implying that

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{2mar}.$$
(24)
This is one-half the angular velocity for the $`1s`$ ground state.

For the $`2p`$ states with $`m=\pm 1`$, cf. Eq. (13), $`S`$ is not identically zero and so we must use the momentum equation (4). We find that

$$𝐩=\pm \frac{\mathrm{}\mathrm{sin}\theta }{2a}\widehat{\varphi }.$$

The coordinates $`r`$ and $`\theta `$ are constant, implying circular orbits about the $`z`$-axis. The orbital angular velocity is given by

$$\frac{d\varphi }{dt}=\pm \frac{\mathrm{}}{2mar}.$$

It will be instructive (for an analysis of Eq. (29) below) to examine the Bohm trajectories resulting from Eq. (10), i.e. ignoring the $`S`$ term. We find that

$$𝐩=\frac{\mathrm{}}{r}\left[\mathrm{sin}\theta (1\frac{r}{2a})+\frac{\mathrm{cos}^2\theta }{\mathrm{sin}\theta }\right]\widehat{\varphi }.$$
(25)
Once again, $`r`$ and $`\theta `$ are constant, implying circular orbits about the $`z`$-axis. However the orbital angular velocity depends upon $`r`$ and $`\theta `$ in a more complicated fashion:

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{mr^2}\left[\frac{r}{2a}1\mathrm{cot}^2\theta \right].$$
(26)
For $`\theta =\pi /2`$, i.e., the $`xy`$ plane, $`\dot{\varphi }=0`$ for $`r=2a`$. This implies that there is a ring of equilibrium points on the $`xy`$ plane at which the electron is stationary. For all other nonzero $`r`$ values on the $`xy`$ plane the electron revolves about the $`z`$-axis: $`\dot{\varphi }<0`$ for $`0<r<2a`$ and $`\dot{\varphi }>0`$ for $`r>2a`$. A ring of stationary points exists on every plane parallel to the $`xy`$ plane: For any fixed $`\theta _0(0,\pi /2)`$, the ring is determined by the relation $`r=2a\mathrm{csc}^2\theta _0`$. The radius of this ring (distance from the $`z`$-axis) is $`r\mathrm{sin}\theta _0=2a/\mathrm{sin}\theta _0`$. The set of all points at which the electron is stationary defines a surface that is generated by revolving the curves

$$z=\pm x\left[\left(\frac{x}{2a}\right)^21\right]^{1/2},x2a,$$
(27)
about the $`z`$-axis. In the region between this surface and the $`z`$-axis, $`\dot{\varphi }<0`$; on the other side of this surface, $`\dot{\varphi }>0`$.

In summary, for each of the hydrogen eigenstates studied above, the Bohm trajectories are circular orbits about the $`z`$-axis, the assumed orientation of the electron spin vector $`𝐬`$. Furthermore, the quantitative behaviour of the angular velocity has been determined in all cases.

We now examine Bohm trajectories for the real hydrogen wavefunctions,

$$\begin{array}{cc}\hfill \psi _{2p_x}& =Nre^{r/2a}\mathrm{sin}\theta \mathrm{cos}\varphi \hfill \\ \hfill \psi _{2p_y}& =Nre^{r/2a}\mathrm{sin}\theta \mathrm{sin}\varphi ,\hfill \end{array}$$
(28)
where $`N=1/\sqrt{32\pi a^5}`$. The probability distributions associated with these wavefunctions are the familiar hydrogen orbitals used in descriptions of organic chemical bonding.

The $`\psi _{2p_x}`$ and $`\psi _{2p_y}`$ wavefunctions are obtained by appropriate linear combinations of the energetically degenerate eigenfunctions $`\psi _{21(\pm 1)}`$ of Eq. (13). Therefore they are also eigenfunctions of the hydrogen atom hamiltonian $`\widehat{H}_0`$ with energy $`E_2`$. From Eq. (9), it follows that $`S=0`$. The cross product of Eq. (10) has components in all three variables, resulting in a system of three coupled ordinary differential equations which must be integrated to find the trajectories.

The system of ODEs associated with the $`\psi _{2p_x}`$ wavefunction is given by

$$\begin{array}{cc}\hfill \frac{dr}{dt}& =\frac{\mathrm{}}{mr}\mathrm{tan}\varphi \hfill \\ \hfill \frac{d\theta }{dt}& =\frac{\mathrm{}}{mr^2}\mathrm{cot}\theta \mathrm{tan}\varphi \hfill \\ \hfill \frac{d\varphi }{dt}& =\frac{\mathrm{}}{mr^2}(1\frac{r}{2a}+\mathrm{cot}^2\theta ).\hfill \end{array}$$
(29)
Note that the $`\varphi `$ DE is identical to Eq. (25). From the first two DEs, we have

$$\frac{dr}{d\theta }=r\mathrm{tan}\theta ,$$
(30)
which is easily integrated to give $`z=r\mathrm{cos}\theta =C`$, in agreement with our earlier analysis. It follows that $`r`$ and $`\theta `$ are constant when $`\varphi =0`$ or $`\pi `$. If $`\theta =\pi /2`$ as well, then $`\dot{\varphi }=0`$ for $`r=2a`$, implying the existence of two equilibrium points at $`(x,y,z)=(2a,0,0)`$ and $`(2a,0,0)`$. At these points, the electron is stationary. In fact, these are two particular cases of an infinity of equilibrium points that are given by the conditions $`\varphi =0,\pi `$ ($`xz`$ plane) and

$$1\frac{r}{2a}+\mathrm{cot}^2\theta =0.$$
(31)
By virtue of the above relation and Eq. (25), the equilibrium points of this system lie on the two curves defined by Eq. (27) in the $`xz`$ plane as well as their reflections about the $`z`$-axis. Each of these points corresponds to the points of highest and lowest “elevation” (from the horizontal $`xy`$ plane) of the familiar dumb-belled level surfaces $`\rho =\psi ^2=C`$ of the orbital. At all of these points, the electron is stationary.

The system of ODEs in (29) may be integrated numerically. However, it is useful to introduce the dimensionless variables $`\xi =r/a`$ and $`\tau =\omega _0t`$, where

$$\omega _0=\frac{E_2E_1}{\mathrm{}}=\frac{3\mathrm{}}{8ma^2}$$
(32)
is the angular frequency associated with the $`n=1`$ to $`n=2`$ transition in the hydrogen atom ($`\omega _01.549\times 10^{16}\text{s}^1`$). The scaled equations become

$$\begin{array}{cc}\hfill \frac{d\xi }{d\tau }& =\frac{8}{3\xi }\mathrm{tan}\varphi \hfill \\ \hfill \frac{d\theta }{d\tau }& =\frac{8}{3\xi ^2}\mathrm{cot}\theta \mathrm{tan}\varphi \hfill \\ \hfill \frac{d\varphi }{d\tau }& =\frac{8}{3\xi ^2}(1\frac{\xi }{2}+\mathrm{cot}^2\theta ).\hfill \end{array}$$
(33)
Figure 1 shows the numerically integrated trajectories for several initial conditions in the $`xz`$ plane. Note that there is a good qualitative agreement between these trajectories and the orbital shapes of the $`2p_x`$ state as depicted in textbook contour plots. (No orbits cross the $`yz`$ plane since it is a nodal surface.) The numerical results confirm that motion is periodic. Angular frequency values are observed to be of the order of $`\omega _0`$. These periodic orbits are stable in the sense of Lyapunov.

The nondimensional and numerical analysis of the $`2p_y`$ case proceeds in a similar fashion. The resulting system of ODEs represents a rotation of the system in Eq. (29) by an angle of $`\pi /2`$ in $`\varphi `$.

To summarize this section, we have shown that the spin-dependent Bohm trajectories for the ground and first excited states of the hydrogen atom are stable periodic orbits. There are no exceptional orbits that deviate from this regularity. For some of the $`2p`$ states, these orbits include families of stationary points that have zero Lebesgue measure in $`𝐑^3`$. These results are intuitively more acceptable than the original Bohmian result that all trajectories associated with a given eigenstate are stationary.

## 3 Trajectories associated with a linear superposition of hydrogenic eigenfunctions

We now examine the Bohm trajectories of an electron with constant spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$ but with a spatial wavefunction $`\psi (𝐱,t)`$ that begins as a linear combination of $`1s`$ and $`2p_0`$ hydrogenic eigenfunctions,

$$\psi (𝐱,0)=c_1\psi _{100}(𝐱)+c_2\psi _{210}(𝐱),$$
(34)
where $`c_1^2+c_2^2=1`$. (The assumption that the electron is in a well-defined “spin up” eigenstate was justified at the end of Section 1.)

The time evolution of this spatial wavefunction, as dictated by the time-dependent Schrödinger equation (8), will be given by

$$\psi (𝐱,t)=c_1\psi _{100}(𝐱)e^{iE_1t/\mathrm{}}+c_2\psi _{210}(𝐱)e^{iE_2t/\mathrm{}}.$$
(35)
This linear combination was chosen in order to examine some of the qualitative features of the trajectories associated with the $`1s`$-$`2p_0`$ transition in hydrogen induced by an oscillating electric field of the form $`\widehat{z}E_0\mathrm{cos}\omega t`$. Indeed, many of the qualitative features of this problem are captured by this model, with the exception of additional oscillations along invariant surfaces due to the oscillating field. We have studied the trajectories of this transition problem in detail and shall report the results elsewhere.

The wavefunction in Eq. (35) is not a linear combination of energetically degenerate states. As such $`S`$, the first term of the momentum in Eq. (4), is nonzero. It can be computed as follows:

$$S=\frac{\mathrm{}}{\psi ^{}\psi }\text{Im}[(\psi )\psi ^{}].$$
(36)
The results are

$`p_r`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{\psi ^{}\psi }}c_1c_2N_1N_2e^{3r/2a}\mathrm{cos}\theta \mathrm{sin}\omega _0t\left(1+{\displaystyle \frac{r}{2a}}\right),`$
$`p_\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{\psi ^{}\psi }}c_1c_2N_1N_2e^{3r/2a}\mathrm{sin}\theta \mathrm{sin}\omega _0t,`$ (37)
$`p_\varphi `$ $`=`$ $`0.`$
Here, $`N_1`$ and $`N_2`$ denote, respectively, the normalization factors of the $`1s`$ and $`2p`$ states and $`\omega _0`$ is defined in Eq. (32). This is the momentum of the original Bohmian formulation $`𝐩=S`$. The angle $`\varphi `$ is constant, implying that there is no orbital motion about the $`z`$-axis.

The second term in the momentum equation (4) contributes only to the $`\varphi `$-momentum:

$$p_\varphi =\frac{\mathrm{}}{\psi ^{}\psi }[Y\mathrm{cos}\theta +X\mathrm{sin}\theta ],$$
(38)
where

$`X`$ $`=`$ $`{\displaystyle \frac{1}{a}}c_1^2N_1^2e^{2r/a}+c_2^2N_2^2r\left(1{\displaystyle \frac{r}{2a}}\right)e^r\mathrm{cos}^2\theta +`$
$`c_1c_2N_1N_2\left(1{\displaystyle \frac{3r}{2a}}e^{3r/2a}\right)\mathrm{cos}\theta \mathrm{cos}\omega _0t,`$
$`Y`$ $`=`$ $`c_2^2N_2^2re^r\mathrm{sin}\theta \mathrm{cos}\theta c_1c_2N_1N_2e^{3r/2a}\mathrm{sin}\theta \mathrm{cos}\omega _0t.`$ (39)
(The terms $`X`$ and $`Y`$ are the $`r`$ and $`\theta `$ components, respectively, of $`\mathrm{log}\psi ^{}\psi `$.) Once again, the spin-dependent momentum term implies orbital motion about the $`z`$-axis.

The $`p_r`$ and $`p_\theta `$ equations of (3) along with $`p_\varphi `$ in (38) yield differential equations in the coordinates $`r`$, $`\theta `$ and $`\varphi `$. Once again, we rewrite these DEs in terms of the dimensionless variables $`\xi =r/a`$ and $`\tau =\omega _0t`$. The net result is the following system:

$`{\displaystyle \frac{d\xi }{d\tau }}`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{3D}}c_1c_2\left(1+{\displaystyle \frac{\xi }{2}}\right)e^{3\xi /2}\mathrm{cos}\theta \mathrm{sin}\tau ,`$
$`{\displaystyle \frac{d\theta }{d\tau }}`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{3D\xi }}c_1c_2e^{3\xi /2}\mathrm{sin}\theta \mathrm{sin}\tau ,`$ (40)
$`{\displaystyle \frac{d\varphi }{d\tau }}`$ $`=`$ $`{\displaystyle \frac{8}{3D\xi }}\left[c_1^2e^{2\xi }+{\displaystyle \frac{1}{64}}c_2^2\xi ^2e^\xi \mathrm{cos}^2\theta +{\displaystyle \frac{3}{8\sqrt{2}}}c_1c_2\xi e^{3\xi /2}\mathrm{cos}\theta \mathrm{cos}\tau \right],`$
where

$$D=c_1^2e^{2\xi }+\frac{1}{32}c_2^2\xi ^2e^\xi \mathrm{cos}^2\theta +\frac{1}{2\sqrt{2}}c_1c_2\xi e^{3\xi /2}\mathrm{cos}\theta \mathrm{cos}\tau .$$
(41)
We first note the following two special cases for the system of DEs in (3) (the primes denote differentiation with respect to $`\tau `$):

1. $`(c_1,c_2)=(1,0)`$ ($`1s`$ state): $`\xi ^{}=\theta ^{}=0`$, and $`\varphi ^{}=8/(3\xi )`$, in agreement with Eq. (7).
2. $`(c_1,c_2)=(0,1)`$ ($`2p_0`$ state): $`\xi ^{}=\theta ^{}=0`$, and $`\varphi ^{}=4/(3\xi )`$, in agreement with Eq. (24).

In these cases, the trajectories are simple circular orbits about the $`z`$-axis, as expected.

More generally, the DEs in $`\xi `$ and $`\theta `$ are decoupled from the $`\varphi `$ DE. From the former two DEs, we have

$$\frac{d\xi }{d\theta }=\xi \left(1+\frac{\xi }{2}\right)\mathrm{cot}\theta .$$
(42)
This separable DE is easily solved to give

$$\xi =\frac{2}{A\mathrm{sin}\theta 1},A=\frac{2+\xi _0}{\xi _0\mathrm{sin}\theta _0}>1,$$
(43)
where $`\xi _0=\xi (0)`$ and $`\theta _0=\theta (0)`$. In a plane with $`z=\xi \mathrm{cos}\theta `$ defining the vertical axis (recall that $`0\theta \pi /2`$), the relation $`\xi (\theta )`$ defines a family of hyperbolae. The asymptotes of the hyperbola in Eq. (43) are given by the rays $`\theta =\theta _1=\text{Sin}^1[\xi _0\mathrm{sin}\theta _0/(\xi _0+2)]`$ and $`\theta =\theta _2=\pi \theta _1`$. If we choose $`\theta _0=\pi /2`$ ($`x`$-axis), then $`\xi _00^+`$, $`\theta _10^+`$, i.e., the hyperbolae flatten as they approach the $`z`$-axis.

Since the $`\xi `$-$`\theta `$ DEs are decoupled from the $`\varphi `$ equation in (3), the electron will remain on the 3D surface obtained by rotating the appropriate $`\xi (\theta )`$ hyperbola in Eq. (43) about the $`z`$-axis. In other words, the hyperboloid of revolution is an invariant set in $`𝐑^3`$ for the electronic trajectory. Plots of three sample trajectories for the case $`c_1=c_2=\frac{1}{\sqrt{2}}`$ are shown in Figure 2. The time interval $`0\tau 50`$ was chosen so that the oscillatory nature of the solutions could be seen.

In Figure 3, the time interval of these solutions has been extended to $`0\tau 1000`$. The invariant surfaces associated with these trajectories can be discerned from the plots.

A more realistic simulation of the oscillating electric field problem is accomplished if we set the coefficients $`c_1`$ and $`c_2`$ in Eq. (34) to be

$$c_1=\mathrm{cos}\omega t,c_2=\mathrm{sin}\omega t,t0,$$
(44)
so that the atom begins in the ground state and proceeds to oscillate between it and the $`2p_0`$ excited state. In this case, the electron trajectories, as determined by by the system of ODEs in (3), will still be constrained to the invariant hyperboloid surfaces of revolution. However, there will be additional oscillatory components in these trajectories due to the periodic behaviour of the $`c_i`$ in Eq. (44). A detailed examination of the possible qualitative behaviour of solutions is beyond the scope of this letter.

## Acknowledgements

We gratefully acknowledge that this research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) in the form of a Postgraduate Scholarship (CC) and an Individual Research Grant (ERV).


text:
# Critical currents at the Bragg glass to vortex glass transition

## Abstract

We present simulations of the transport properties of superconductors at the transition from the Bragg glass (BG) to the vortex glass (VG) phase. We study the frustrated anisotropic 3D XY model with point disorder, which has been shown to have a first order transition as a function of the intensity of disorder. We add an external current to the model and we obtain current-voltage curves as a function of disorder at a low temperature. We find that the in-plane critical current has a steep increase at the BG-VG transition, while the $`c`$-axis critical current has a discontinous jump down, this later result in agreement with the first-order character of the transition.

The study of the vortex phase diagram of anisotropic superconductors with point disorder has been of interest since the discovery of high $`T_c`$ superconductivity. It is now clear that at low temperatures and low magnetic fields there is a “Bragg glass” phase (BG) with an elastically distorted vortex lattice giamarchi . This vortex lattice undergoes a first order melting transition to a vortex liquid (VL) when increasing the temperature $`T`$. At higher magnetic fields and low $`T`$ there is a disordered vortex state, the “vortex glass” (VG) vglass ; bokil ; kosterlitz ; olsson . A disorder driven transition from the BG to the VG has been proposed giamarchi  to occur when increasing the magnetic field $`H`$ or the disorder. Experimental observations with increasing $`H`$ at low $`T`$, that have been attributed to the BG-VG transition, are: the destruction of the Bragg peaks in neutron scattering neutron , a dip in the differential resistance safar , the onset of the “second magnetization peak” khay ; nishizaki , or a jump in the Josephson plasma resonance gaifu . Numerical simulations in particle-like models with random pinning have found a transition from an ordered lattice to a disordered lattice ryu ; zimanyi  when increasing particle density (i.e., the magnetic field). Recently, Monte Carlo simulations in the frustrated 3D XY model with disorder have stablished that the BG-VG transition is of first-order type hu ; teitel .

Several of the experimental evidences of the BG-VG transition are directly or indirectly related to transport properties. For example, the most common determination of the BG-VG line is through the onset of a “second magnetization peak” khay ; nishizaki , which is attributed to an steep increase of the in-plane critical currents. From the point of view of transport, the BG has zero linear resistivity ($`\rho _{\mathrm{lin}}=0`$) giamarchi ; kosterlitz . The VG phase was originally proposed to have zero resistivity vglass . However, studies of the XY gauge glass model with finite screening ($`\kappa <\mathrm{}`$) have found that VG is unstable in $`d=3`$ bokil ; kosterlitz , and the disordered phase is a frozen vortex liquid that has a very small finite resistivity. Only for the $`\kappa \mathrm{}`$ disordered XY and gauge glass models the VG phase is stable in $`d=3`$ at low $`T`$, in which case $`\rho _{\mathrm{lin}}=0`$ olsson . Therefore, it is important to understand how the the non-linear current-voltage (IV) curves and their critical currents change across a disorder driven first-order transition. We present here numerical calculations of both the in-plane and the $`c`$-axis IV curves across the disorder driven BG-VG transition in the $`\kappa =\mathrm{}`$ frustrated XY model used in hu ; teitel  (which is the only model where the first-order nature of the transition has been clearly determined so far). We find that the $`c`$-axis critical current has a discontinuous jump at the BG-VG transition, which could be tested experimentally.

The hamiltonian of the frustrated 3D XY model is:

$$H=\underset{𝐫}{}\underset{\mu =\widehat{𝐱},\widehat{𝐲}}{}J_{𝐫\mu }V(\theta _\mu (𝐫))\frac{J}{\mathrm{\Gamma }}\mathrm{cos}(\theta _{\widehat{𝐳}}(𝐫))$$

where $`\theta _\mu (𝐫)=\theta (𝐫+\mu )\theta (𝐫)A_\mu (𝐫)`$ and $`A_\mu (𝐫)=\frac{2\pi }{\mathrm{\Phi }_0}_{𝐫a}^{(𝐫+\mu )a}𝐀𝑑𝐥`$ and $`𝐫=(n_x,n_y,n_z)`$ defines a lattice site in a cubic grid. To minimize the pinning of the numerical grid, we have chosen $`V(\theta _\mu (𝐫))=R_0R_1\mathrm{cos}(\theta _\mu (𝐫))R_2\mathrm{cos}(2\theta _\mu (𝐫))`$ and adjusted the coefficients $`R_0,R_1`$ and $`R_2`$ as in koshelev . We model uncorrelated random point pinning in the xy plane hu ; teitel  with $`J_{𝐫\mu }=J(1+pϵ_{r,\mu })`$, $`\mu =\widehat{𝐱},\widehat{𝐲}`$, where $`ϵ_{r,\mu }`$ are independent random variables with $`ϵ_{r,\mu }=0`$ and $`ϵ_{r,\mu }^2=1`$. In the presence of an external magnetic field $`H`$ along the $`\widehat{z}`$ direction, we have $`A_x(𝐫)A_x(𝐫+𝐲)+A_y(𝐫+𝐱)A_y(𝐫)=2\pi f`$, and we take $`f=Ha^2/\mathrm{\Phi }_0=\frac{1}{24}`$. We model the current between two grid points with the RSJ model dd :

$$I_\mu (𝐫,t)=S_\mu (𝐫,t)+N_\mu (𝐫,t)+\eta _\mu (𝐫,t)$$

with the superconducting current $`S_\mu (𝐫,t)=\frac{2e}{\mathrm{}}\frac{H}{\theta _\mu (𝐫)}`$, the normal current $`N_\mu (𝐫,t)=\frac{\mathrm{}}{2eR_\mu }\frac{\theta _\mu (𝐫)}{t}`$ with $`R_\mu =R_{ab}`$ for $`\mu =\widehat{𝐱},\widehat{𝐲}`$ and $`R_\mu =\mathrm{\Gamma }R_{ab}`$ for $`\mu =\widehat{𝐳}`$, and the thermal noise fluctuations have correlations $`\eta _\mu (𝐫,t)\eta _\mu ^{}(𝐫^{},t^{})=(2k_BT/R_\mu )\delta _{\widehat{\mu },\widehat{\mu }^{}}\delta _{𝐫,𝐫^{}}\delta (tt^{})`$. We take periodic boundary conditions in all directions with a fluctuating twist $`\alpha _\mu `$ such that $`A_\mu (𝐫,t)=A_\mu ^0(𝐫)\alpha _\mu (t)`$. This allows to obtain the voltage in each direction as $`V_\mu =\frac{\mathrm{}}{2e}\frac{d\alpha _\mu (t)}{dt}`$ periodic . We consider current conservation in each node: $`_\mu I_\mu (𝐫)I_\mu (𝐫\mu )=0,`$ and we fix the total current in each direction, consistently with the periodic boundary condition, by periodic :

$$I_\mu ^{\mathrm{ext}}=\frac{\mathrm{}}{2eR_\mu }\frac{d\alpha _\mu }{dt}+\frac{1}{L_xL_yL_z}\underset{𝐫}{}\left[S_\mu (𝐫)+\eta _\mu (𝐫,t)\right]$$

We integrate with a second order Runge-Kutta stochastic algorithm with $`\mathrm{\Delta }t=0.1\mathrm{}^2/(4e^2R_{ab}J)`$ and we average over an interval of $`10^510^6\mathrm{\Delta }t`$. We consider system sizes $`48\times 48\times L_z`$ with $`L_z=1232`$ and anisotropy $`\mathrm{\Gamma }=40`$. Most of the results shown are for $`L_z=32`$ and for a single disorder realization. We have considered other realizations of disorder for some values of $`T`$ and $`p`$ and have found similar behavior.

We first study the model at equilibrium (in the absence of external currents). At low disorder, $`p=0.01`$, we find a transition upon increasing $`T`$ from a vortex lattice to a vortex liquid. We have calculated the structure factor $`S(𝐤)=(1/Nf^2)|_{𝐫_{},z}n_z(𝐫_{},z)e^{i𝐤.𝐫_{}}|^2`$, where $`n_z`$ is the vorticity in the $`xy`$ planes and $`N=L_xL_yL_z`$. In Fig. 1(a) we plot the intensity of the structure factor, $`S(𝐆)`$, at one of the first reciprocal lattice vectors $`𝐆`$, showing the melting of the vortex lattice at $`T_c0.28`$. The resistivity is calculated by driving the system with a very small current and obtaining voltage. We show in Fig.1(a) the $`c`$-axis resistivity $`\rho _c`$ vs. $`T`$ which drops sharply at $`T_c`$; we find similar behavior in $`\rho _{ab}`$ (not shown). At a low temperature $`T=0.1J/k_B`$, we vary disorder and study the behavior of the lattice order parameter, $`S(𝐆)`$, which is shown in Fig. 1(b) for two system sizes. We find that $`S(𝐆)`$ drops at $`p_c0.11`$, indicating a transition from the BG to the VG. Typical structure factors $`S(𝐤)`$ for $`p<p_c`$ and $`p>p_c`$ are shown in the insets. Repeating this analysis for other $`T`$, we have found that $`p_c(T)`$ is approximately independent of $`T`$ for low temperatures. Similar phase boundary $`p_c(T)`$ for the BG phase were reported in hu ; teitel . For very low temperatures, both hu  and teitel  agree that there is a first-order transition from the BG phase to the VG phase in this model. (At higher $`T`$ there is disagreement on the existence of a “vortex slush” phase).

From now on, we focus on the disorder driven transition at low $`T`$, where hu ; teitel  agree on the phase diagram. We fix a low temperature $`T=0.1J/k_B`$ and we vary the disorder strength $`p`$. In Fig. 2(a) we show the calculated in-plane IV curves for different $`p`$.

They are obtained by applying a uniform current along the $`x`$ direction, $`I_\mu ^{\mathrm{ext}}=I_{ab}\delta _{\mu ,x}`$, and calculating the average voltage along $`x`$, $`V_{ab}=d\alpha _x/dt`$. In Fig. 2(b) we show the $`c`$-axis IV curves. In this case we apply a uniform current along the $`\widehat{z}`$ direction (parallel to the magnetic field), $`I_\mu ^{\mathrm{ext}}=I_z\delta _{\mu ,z}`$, and we calculate the voltage $`V_z=d\alpha _z/dt`$. Before doing further analysis it is important to discuss the effect of a small driving current in the BG-VG transition. We find that currents much smaller than the critical currents have a negligible effect on $`p_c`$. In Fig. 3 we show $`S(𝐆)`$ vs. $`p`$ when the driving current in a given direction is of the order of the corresponding critical current. In the case of a $`z`$-direction current $`I_z=0.010I_c^z`$, there is a very small effect with a critical $`p_c`$ slightly smaller than the $`I=0`$ case. In the case of an in-plane current $`I_{ab}=0.002I_c^{ab}`$, we find that there is a first drop of $`S(𝐆)`$ to a finite value at a $`p<p_c`$. This corresponds to the situation when the applied current $`I_{ab}`$ equals the critical current for that $`p`$. Then a second drop of $`S(𝐆)`$ to the very small values of the VG phase occurs near the critical $`p_c`$ but also at a slighlty smaller value. In both cases, one can conclude that finite currents of the order of the critical current have a small effect on the value of $`p_c`$, lowering it in about $`\mathrm{\Delta }p_c0.005`$.

Let us now discuss the in-plane transport at the transition. In Fig.2(a) we show with open symbols the IV curves for $`p<p_c`$ and with filled symbols the IV curves for $`p>p_c`$ (where $`p_c=0.11`$ is the zero current critical disorder). We observe that in both phases the IV curves are strongly nonlinear and that they tend to zero resistivity for low currents ($`V/I0`$). In order to observe if there is a change when crossing the BG-VG line, we fix a value of the current $`I_{ab}`$, slightly above the critical currents, and we calculate the voltage $`V_{ab}`$ as a function of disorder strength $`p`$. This is shown in Fig.4(a). We see that at $`p_c`$ there is a jump down in $`V_{ab}`$.

This jump is observed for different values of $`I_{ab}`$. The location of the jump does not depend strongly on $`I_{ab}`$. The experiment of safar  finds a clear decrease of the differential resistance, measured for a current above the critical current, when increasing magnetic field at low temperatures in YBaCuO samples. This is in agreement with our result of Fig. 4(a) for the behavior of the in-plane voltage at the BG-VG transition. One can interpret this jump down in voltage as a consequence of a jump up in the critical current. It is difficult to define a critical current, since at $`T>0`$ there is always a non-zero voltage for any current. Experimentally, the “critical current”, $`I_c`$, is obtained from the IV curves as $`V(I_c)=V_{th}`$, with the threshold voltage $`V_{th}`$ small enough. If $`I_c`$ is almost independent of $`V_{th}`$, then it corresponds to a crossover current separating two regimes in the IV curve. In our case, $`I_c`$ separates the ohmic flux flow regime from the low current nonlinear regime which has $`\rho _{\mathrm{lin}}(I0)=0`$ for the two phases. (In the case of the VG with finite screening, $`\kappa <\mathrm{}`$, $`I_c`$ would separate two ohmic regimes). To obtain the $`ab`$ critical current, $`I_{ab}^c`$, we choose different values of $`V_{th}`$ slightly above our noise value of voltage. One of the $`V_{th}`$ used is shown in Fig.2(a) with a dashed line. The obtained critical currents are plotted in Fig.5(a) as a function of $`p`$ (we see that different $`V_{th}`$ give smilar values of $`I_c^{ab}`$, within the error bars). We observe in the plot that there is a steep increase of $`I_c^{ab}`$ at $`p_c`$.

In the magnetization experiments in YBaCuO nishizaki  the BG-VG transition line has been associated to the onset of the second magnetization peak, which has been related to an steep increase of the in-plane critical currents. Our result in Fig 5(a) is in good agreement with this interpretation, particularly because the value of $`\mathrm{\Gamma }`$ we use corresponds to YBaCuO. The more anisotropic BiSrCaCuO samples show experimentally a much stronger increase in the second magnetization peak khay . In this very anisotropic case the increase in critical current has also been interpreted as a decoupling transition colson .

The 3D XY model is very adequate to the description of the $`c`$-axis transport, since the Josephson coupling between planes is treated exactly. The $`c`$-axis current-voltage charactestics, $`V_z`$ vs. $`I_z`$, are shown in Fig.2(b) for different values of $`p`$ at $`T=0.1J/k_B`$. The IV curves for $`p<p_c`$ are plotted with open symbols and the IV curves for $`p>p_c`$ are plotted with filled symbols. We observe that in all cases the IV curves show a well-defined critical current where the voltage drops steeply. The IV curves can be clearly separated in two different sets of curves, one set for $`p<p_c`$ with high critical currents and another set for $`p>p_c`$ with low critical currents. This can also be observed if for a fixed current $`I_z`$ we vary $`p`$ and we calculate $`V_z`$, as we show in Fig.4(b). For very low $`I_z`$ we find $`V_z0`$ for all $`p`$, consistent with both phases being superconducting. For intermediate values of $`I_z`$ we find that $`V_z0`$ for $`p<p_c`$ with a sharp jump to finite voltage $`V_z`$ at $`p_c`$. For values of $`I_z`$ above the critical currents, there is a finite voltage $`V_z`$ in both phases, with a clear discontinuous jump up at $`p_c`$. We obtain the critical currents $`I_c^z`$ with a voltage criterion $`V_{th}`$, which is shown in Fig.2(b) with a dashed line. It is obvious in this case that the values obtained for $`I_c^z`$ are insesitive to the choice of $`V_{th}`$. The $`c`$-axis critical currents are plotted as a function of disorder strength $`p`$ in Fig. 5(b). We observe in the plot that $`I_c^z`$ has a clear discontinuous jump down at $`p_c`$. This result is consistent with the Josephson plasma resonance measurements of gaifu , which are sensitive to the $`c`$-axis Josepshon coupling, where a sharp decrease of the plasma resonance frequency was observed at the transition. Recently, the $`c`$-axis Josephson critical current has been measured experimentally in BiSrCaCuO samples varying magnetic field at the melting transition line, the BG-VL line ooi . A jump down of the critical current was found, although at a field slightly lower than the equilibrium case. It will be interesting if a similar experiment could be carried out across the BG-VG line, since our results imply that $`c`$-axis transport would give a better evidence (compared to the $`ab`$-plane transport) of the first-order nature of this transition.

We acknowledge financial support from Conicet, CNEA, ANPCYT (PICT99-03-06343) and Fundación Antorchas (Proy. 14116-147). A.D.H. also ackowledges support from the Centro Latinoamericano de Física.


text:
# Methods of calculation of a friction coefficient: Application to the nanotubes

## Abstract

In this work we develop theoretical and numerical methods of calculation of a dynamic friction coefficient. The theoretical method is based on an adiabatic approximation which allows us to express the dynamic friction coefficient in terms of the time integral of the autocorrelation function of the force between both sliding objects. The motion of the objects and the autocorrelation function can be numerically calculated by molecular-dynamics simulations. We have successfully applied these methods to the evaluation of the dynamic friction coefficient of the relative motion of two concentric carbon nanotubes. The dynamic friction coefficient is shown to increase with the temperature.

Thanks to recent developments in nanotechnology, the hope is high to build mechanical devices on the scale of the nanometer. For this purpose, it is important to determine the mechanical properties and, especially, the friction forces in such nanodevices. In this work, we are going to focus on systems of carbon nanotubes. Our particular interest is for a system that Cumings and Zettl Zettl  observed experimentally with a TEM. They fixed an edge of a multiwalled nanotube to a surface and opened the other edge, they then extracted inner layers from the core for several nanometers and released them. They observed a full retraction of the inner layers and furthermore, they could conclude that the multiwalled nanotubes are self cleaning since amorphous carbon due to the opening of the edge are inside the tube and also do not have any wear, structural change or fatigue after several extraction and retraction processes. This experiment suggest us that the multiwalled nanotubes can be promising systems for future nanometric mechanical parts such as springs, gears or even motors.

A short time after the work of Cumings and Zettl, Zheng and Jiang Zheng  estimated the frequency of the oscillations in this system to be of the order of GHz. This result is also very interesting since in the macroscopic world moving parts with such frequencies does not exist at the present time. But as in the macroscopic world moving parts have friction forces which hinder the motion and dissipate energy. We hence have to know the importance of these forces before conceiving such devices. Our work will thus focus on, firstly methods of calculation of the friction and secondly an application of these theories to the multiwalled nanotubes.

The plan of the paper will be as follows, we are first going to introduce the theoretical framework for the description of the mechanics of nanotubes. We will then solve the classical equations of motion and calculate the dynamic friction coefficient by molecular-dynamics simulation and by the autocorrelation-function method developed by Jarzynski, Berry and Robbins Jarzynski ; Berry .

We have depicted in Fig. 1 two sliding nanotubes. $`R`$ is the distance between the centers of mass of each nanotube.

The Hamiltonian of the system of two nanotubes is given by

$$H=T_1+T_2+V_{\mathrm{TB}}^{(1)}+V_{\mathrm{TB}}^{(2)}+\frac{1}{2}\underset{ij}{}V_{\mathrm{LJ}}\left(𝐫_i^{(1)}𝐫_j^{(2)}\right)$$
(1)
where $`T_1`$ and $`T_2`$ are respectively the total kinetic energies of the inner and outer tubes, $`V_{\mathrm{TB}}^{(1)}`$ and $`V_{\mathrm{TB}}^{(2)}`$ are the total Tersoff-Brenner potentials of both tubes. We use the set of parameters which optimizes the geometric structure Brenner . The last term is the interaction potential between the tubes, which is taken as a 6-12 Lennard-Jones potential calculated by Lu and Yang Lu  and widely used for molecular-dynamics simulations of nanotubes Zheng ; Lu ; Lu2 ; Legoas .

The calculation of the friction coefficient will be done by using an adiabatic approximation since there are two different time scales in the system. The fast time scale is the one of the vibrations of the atoms which is of the order of the $`fs`$ and the slow time scale is the one of the motion of the mass centers of the tubes which is of the order of the $`ps`$ in our case. Thanks to this difference in time scales, we can use a method developed by Jarzynski, Berry and Robbins Jarzynski ; Berry  which is based on the fact that the following quantity is an adiabatic invariant Ott ; Goldstein ,

$$\mathrm{\Omega }(E,ϵt)=𝑑𝐳\theta \left[EH(𝐳,ϵt)\right]$$
(2)
where $`ϵ`$ is a small parameter, E the energy of the fast system and $`H`$ is an ergodic adiabatic Hamiltonian. This means that $`H`$ slowly evolves with time and that, at any instantaneous time, $`H`$ produces trajectories which explore their energy shells ergodically and chaotically. $`𝐳=(𝐪,𝐩)`$ represents the fast phase-space variables. The smallness of the parameter $`ϵ`$ guarantees that the motion of vibration is fast compared to the motion of the mass centers. The other important assumption is that the force on the slow system can be expressed as an average over the fast degrees of freedom since the time scale of the fast system is much smaller. We hence write the force as

$$𝐅=𝑑𝐳\rho H$$
(3)
where $`\rho `$ is the probability density of the fast system and $`=/𝐑`$ is the derivative with respect to the slow variable. In our specific case $`𝐑`$ is the distance between the mass centers of the tubes. Using the fact that there are two different time scales, $`\rho `$ is expanded as

$$\rho (𝐳,t)=\rho _0(𝐳,𝐑)+ϵ\rho _1(𝐳,𝐑,t)+\mathrm{}$$
(4)
where the first order term $`\rho _0`$ depends only on the slow time which is regarded as the parameter $`𝐑`$. The evolution of $`\rho `$ is given by the Liouville equation. Injecting the expansion (4) in the Liouville equation gives a hierarchy of equations which can be solved using a microcanonical distribution for $`\rho _0`$. Having calculated $`\rho `$ to the first order in $`ϵ`$ one can find the following reaction force,

$$𝐅=E𝚪\dot{𝐑}$$
(5)
where the first order term is the Born-Oppenheimer force with the energy of the fast system $`E`$ and the second term the friction force where the velocity $`\dot{𝐑}`$ of order $`ϵ`$ appears. The dynamic friction coefficients $`𝚪`$ are given by

$$\mathrm{\Gamma }_{ij}=\frac{1}{_E\mathrm{\Omega }}_E\left(_E\mathrm{\Omega }_0^{\mathrm{}}𝑑\tau C_{ij}(\tau )\right)$$
(6)
where $`_E`$ is the partial derivative with respect to the energy of the degrees of freedom of vibration $`E`$ and $`\mathrm{\Omega }(E,𝐑)`$ is defined by Eq. (2). The autocorrelation function is given by the following expression

$$C_{ij}(\tau )=_i\stackrel{~}{H}_\tau _j\stackrel{~}{H}_{E,𝐑}$$
(7)
where $`\stackrel{~}{H}=H(z,𝐑)E(𝐑)`$ represents the fluctuations of the Hamiltonian. We can simplify the expression of the friction in Eq. (6) by using the Boltzmann equation $`S=k_\mathrm{B}\mathrm{ln}_E\mathrm{\Omega }`$ and the thermodynamic relation $`E/S=T`$. Moreover, the adiabatic parameter $`𝐑`$ is one-dimensional in our system since the relative motion of the tubes is mainly in the longitudinal direction. We hence get

$$\mathrm{\Gamma }=\left(\beta +_E\right)_0^{\mathrm{}}𝑑\tau C(\tau )\beta _0^{\mathrm{}}𝑑\tau C(\tau )$$
(8)
with $`\beta =1/k_\mathrm{B}T`$. In Eq. (8), we have used the result that the second term with $`_E`$ is a negligible correction of the order of the inverse of the system size. The autocorrelation function is given by

$$C(\tau )=F_{\mathrm{LJ}}(\tau )F_{\mathrm{LJ}}(0)_{E,R}F_{\mathrm{LJ}}_{E,R}^2$$
(9)
where $`F_{\mathrm{LJ}}`$ is the total Lennard-Jones force in the transverse direction,

$$F_{\mathrm{LJ}}=\frac{1}{2}\underset{ij}{}\frac{V_{\mathrm{LJ}}}{R}\left(𝐫_i^{(1)}𝐫_j^{(2)}\right)$$
(10)
We observe that Eq. (8) satisfies the fluctuation-dissipation theorem relating the fluctuating force acting on a system in a bath to the friction kernel in the generalized Langevin equation:

$$k_\mathrm{B}T\zeta (t)=f(t)f(0)$$
(11)
with $`f=F_{\mathrm{LJ}}F_{\mathrm{LJ}}`$. Eq. (8) is hence obtained in the Markovian limit where the friction kernel (11) decays on a short time scale such that $`\zeta (t)=\zeta _0\delta (t)`$.

We are now going to show the results of the molecular-dynamics simulation. We focus our numerical studies on a specific nanotube namely (5,0)@(15,0) which is commensurate in the classification of Ref. Crespi . The inner tube has $`N_1=60`$ atoms and a length of $`l_1=10.9\text{Å}`$. The outer tube has $`N_2=240`$ atoms and a length of $`l_2=15.3\text{Å}`$. The relative mass is hence $`\mu =576`$ amu. The molecular-dynamics simulation is straightforward. We first extract the internal nanotube such as the initial distance between both mass centers is $`R_0=6.0\text{Å}`$. The tubes centers of masses have vanishing initial velocity and angular momentum. We then start the simulation and solve the Newtonian dynamics with a fourth-order Runge-Kutta integration scheme in the microcanonical ensemble. The initial temperature of the system is $`300K`$, which is introduced by a Maxwell-Boltzmann distribution of the initial velocities and a short early time relaxation. We use a time step of $`0.1`$fs.

First, we can predict the energy to be dissipated by the sliding motion and thus the increase of the temperature of the system. Calculating the difference of energies between the minimum of the total Lennard-Jones potential $`V_{\mathrm{LJ}}(R)`$ and the initial value $`V_{\mathrm{LJ}}(R_0)`$ gives us $`\mathrm{\Delta }E_s=0.623eV`$. Using $`\mathrm{\Delta }E=3Nk_BT`$ where $`N=N_1+N_2`$ we hence should have an increase in temperature of $`\mathrm{\Delta }T=8.04K`$.

We observe in Fig. 2 a double exponential decay in the motion. The origin of the double exponential decay holds in the fact that, in the initial configuration, a part of the inner tube is outside the outer tube. In this case, the friction is higher because the inner tube interacts with the edge of the outer tube and bounces during the reentry. It is only after a while that the inner tube oscillates inside the outer tube. We see that the oscillations are damped as expected. We can use a simple model to describe the dynamics when the inner tube remains inside, in which case the potential energy of the relative motion between both tubes can be approximated by a harmonic potential. We hence have a Newton equation of the form,

$$\mu \ddot{R}=\mathrm{\Gamma }\dot{R}kR$$
(12)
which can be solved to show that the successive maxima in position $`R`$ and energy $`E=\frac{1}{2}\mu \dot{R}^2+\frac{1}{2}kR^2`$ decay as

$$|R|_{\mathrm{max}}e^{\mathrm{\Gamma }t/2\mu }\text{and}E_{\mathrm{max}}e^{\mathrm{\Gamma }t/\mu }$$
(13)
We hence expect to have a ratio of 2 between the damping rates of the maxima in energy and position. By performing fits of the curves in Fig. 2, we get respectively $`0.00072`$/ps and $`0.0015`$/ps, in good agreement with the expectation. Multiplying by the relative mass $`\mu `$ we finally get the friction coefficient $`\mathrm{\Gamma }=0.84\pm 0.02`$ amu/ps.

The next step is to calculate the friction coefficient with the autocorrelation-function method. According to Eqs. (8)-(10), the friction coefficient is given by the time integral of the autocorrelation function of the force acting on the slow system with the distance between the centers of mass fixed as a parameter. The calculation of the autocorrelation function is also performed by molecular-dynamics simulation except that the distance between the centers of mass is kept constant by compensating at each time step the motion of each atom by the motion of the mass center of their respective nanotube. We suppose that our system is ergodic and mixing (meaning that its autocorrelation functions vanish asymptotically at long times). Since we have nearly one thousand degrees of freedom this should be a good approximation. We can thus use the ergodic theorem and replace the ensemble average in Eq. (9) by a time average.

Figure 3 depicts the results of the numerical calculation of the autocorrelation-function and its time integral giving asymptotically the friction coefficient. With this autocorrelation-function method, we find the value $`\mathrm{\Gamma }=0.82`$ amu/ps for the friction coefficient, in excellent agreement with the value obtained by direct molecular-dynamics simulation. The great advantage of the autocorrelation function method is that the calculation is now about 50 times faster. Furthermore the autocorrelation function method is much more precise since it does not involve fitting to get the friction coefficient.

Furthermore, the autocorrelation-function method allows us to obtain a theoretical formula for the friction coefficient. We first notice that the autocorrelation function decays on the time scale of the vibrations of the atoms in the carbon nanotubes. This time scale is of the order of the inverse of the Debye frequency $`\omega _\mathrm{D}`$ so that the correlation time should be given approximately by $`t_\mathrm{C}2\pi /\omega _\mathrm{D}`$. Since the Debye temperature of the nanotube is $`\mathrm{\Theta }_\mathrm{D}1000`$K Debye , we hence expect the autocorrelation function to decay to zero over a correlation time approximately equal to 50fs, which is indeed confirmed by Fig. 3. Moreover, the autocorrelation function at time zero is equal to the variance $`\mathrm{\Delta }F^2`$ of the force. We can therefore approximate the time integral of the autocorrelation function and obtained the following approximate formula for the friction coefficient:

$$\mathrm{\Gamma }\frac{1}{k_\mathrm{B}T}\frac{1}{2}\frac{2\pi }{\omega _\mathrm{D}}\mathrm{\Delta }F^2$$
(14)
Thanks to this formula, the evaluation of the friction coefficient reduces to the evaluation of the variance of the fluctuating force. On this ground, we now investigate the dependence of the friction coefficient on temperature. We have depicted in Fig. 4 the result of the numerical calculation of the standard deviation of the force.

We see that the standard deviation of the force increases linearly with the temperature, hence the friction coefficient also increases linearly with temperature: $`\mathrm{\Gamma }T`$. The fact that friction increases with temperature in our system can be understand by the fact that the dissipation of energy of the slow motion is due to the coupling with the vibrating degrees of freedom. Increasing the temperature results into larger vibrational motion and thus into higher friction. Furthermore, the calculation shows that edge effects give to the small system we here considered, a larger friction than to a system of double length. Then the friction increases in systems from double to longer lengths, the coefficient remaining of the same order of magnitude as in the small system here studied.

In this work, we have developed theoretical and numerical methods to calculate the dynamic friction coefficient in the one-dimensional sliding motion of two concentric carbon nanotubes. We focused our numerical studies on a commensurate double-walled nanotube. We obtained by a molecular-dynamics simulation the energy damping and thus the increase in temperature of the system. The simulation shows the importance of edge effects at early times: a significantly larger damping is due to bouncing effects when the inner nanotube reenters. This effect disappears at long times when the inner tube moves inside the outer tube. The long time behavior is in good agreement with the one-dimensional model we proposed. We have then calculated the friction with the autocorrelation-function method. We observed that the autocorrelation function decreases to zero with a time of the order of 50fs which is a value in agreement with the experimental value of the Debye temperature. The integral of the autocorrelation function gives the same result as the direct molecular-dynamics simulation with much less calculation time, actually a factor 50, and with better precision since we avoid the use of fitting methods. We can thus conclude that the autocorrelation-function method is very efficient for the calculation of friction coefficients. Moreover, we have obtained a formula which relates the dynamic friction coefficient to the Debye frequency and the standard deviation of the fluctuations of the forces between the nanotubes and determined that the dynamical friction coefficient increases with the temperature.

The Authors thank Pablo Jensen for helpful discussions. This research is financially supported by the National Fund for Scientific Research (FNRS Belgium).