text:
# 1 a) German mortality functions and Gompertz law, for 1875, 1905, 1933, and 1996. b) Remaining life expectancy for German men at age 65, and simulation results from the simple model presented later.

Simple tools for extrapolations of human mortality in rich countries

Dietrich Stauffer

Institute for Theoretical Physics, Cologne University, D-50923 Köln, Euroland

Abstract: Suitable assumptions for the Gompertz mortality law take into account the break in the time development observed recently by Wilmoth et al. They show how a drastic reduction in the birth rate and improved living conditions lead to a drastic increase in the fraction of old people in the population, and how immigration of half a percent of the population per year can mostly stop this increase.

Human life expectancies increased and birth rates decreased, in the 20th century in most rich countries, and this may happen soon elsewhere . Predictions for the future \[2-4\] are usually extrapolations of the trends of the past and thus are problematic if extended over centuries as in ref.4. Recent analysis of the age of the oldest people in a finite population versus time , of survival probabilities as a function of life expectancy at birth , or of survival probabilities as a function of survival probability up to age 40  raise the possibility of sharp breaks in these curves. By definition, no extrapolation beyond the next (unknown) break is possible, if these breaks are infinitely sharp (if infinitely accurate data would be available). Now we try to describe these effects by a simple model with a limited number of parameters, and show the versatility of this model by an extrapolation of the fraction of old people under various hypotheses. While some aspects are motivated by German data, the one-page computer program could easily be adapted to other rich countries with birth rates below replacement level and is available from stauffer@thp.uni-koeln.de.

Human adult mortalities between the ages $`x`$ of 30 and 90 years follow roughly an exponential increase $`\mathrm{exp}(bx)`$, known as the Gompertz law since the 19th century. Much of the trends of the past centuries in various countries were summarized by the universality law :

$$\mu /b=Ae^{b(xX)}$$

where $`\mu `$ is the mortality function, $`A`$ and $`X103`$ years are universal over the whole human species, and only the slope $`b`$ increases with increasing progress. The time unit is one year throughout. (Also Thatcher  noted that human mortality at about 100 years is roughly constant.) We take $`b`$ to increase linearly with time from 0.07 in 1821 to 0.093 in 1971, and from then on to remain constant. Fig.1a compares the mortality function of German women with this simple assumption. We see that our analysis is valid only at older ages and only for the last decades of the 20th century in rich countries when child mortality was already low. Fig.1b compares the life expectancy at age 65 for German men with simulations of the model explained below.

(In the above equation, the mortality function $`\mu =d\mathrm{ln}S(x)/dx`$ is also called the hazard factor or force of mortality  and is approximated by $`\mathrm{ln}[S(x1/2)/S(x+1/2)]`$ where $`S(x)`$ (often denoted as $`l(x)`$) is the probability to survive from birth to age $`x`$. For centenarians, downward deviations of $`\mu `$ from the Gompertz law have been reported  but an extrapolation to zero systematic errors (from age overstating) seems lacking.)

However, these assumptions are still wrong since they lead for increasing $`b`$ (increasing progress) to a decrease of the average maximum age in a finite population. (This maximum age $`x_{max}`$ is defined as the age at which the reciprocal survival probability $`1/S(x_{max})`$ reaches the population size ; we ignore the growth of the population with time and also the variation of the mortality parameters during the lifetime of an individual.) In reality, $`x_{max}`$ increased with increasing progress during the last century , page 46 in .

To allow for $`x_{max}`$ to increase with calendar year $`t`$, we allow the characteristic age $`X`$ to increase with time after 1971, by generalizing $`X=103`$ to $`X=103+0.15(t1971)`$. This generalization is already included in Fig.1 and slightly separates the lowest straight line from the three higher (earlier) ones. Now $`x_{max}`$ increases with increasing $`t`$, similar to reality (page 46 in .)

The two effects, the increase of $`b`$ before the year 1971, and the increase of $`x`$ thereafter, describe the two empirical observations reviewed recently by Yashin et al : During the first half the the 20th century the survival curves $`S(x)`$ became more rectangular (“compression of mortality”), but later this rectangularization stopped and is replaced by a shifting of the old-age survival curve towards older ages at constant shape. The interquartile death rate is the number of years between $`S(x)=0.75`$ and $`S(x)=0.25`$. In my calculations it first declines and then stays constant similar to Swedish reality (Fig.4a of Wilmoth and Horiuchi , if deaths before age 30 are ignored).

For the birth rate we assume that pregnancies happen at ages 20 to 40 only, with equal probability in each of these 20 years. We first adjust a birth rate $`m=2.17`$ such that for $`b=0.07`$ (a realistic value for earlier centuries ) a stationary population is reached. (Note that “births” here count only babies reaching adult age; child mortality was ignored here from the beginning.) The birth rate then increased linearly in time to 2.2 for the year 1971. Then very rapidly it decreases to 1.4 and stays there, according to a hyperbolic tangent formula

$$m(t)=1.80.4\mathrm{tanh}[(t1971)/3]$$

compared in Fig.2 with West German reality.

Thus, compatible with Wilmoth et al , the year 1971 is taken as both the year when rectangularization was replaced by shifting in $`S(x)`$, and then the birth rate had its strongest decline. These sharp changes in the higher time derivatives may be missing in traditional demographic extrapolations \[2-4\] but could be easily implemented in a Fortran program of about 50 lines. Thus these methods could be used in more realistic fits and extrapolations for specific countries.

Figure 3 illustrates from these assumptions the future development for the fraction of “old” people, the limit taken at 65. About this “age quake” much has been talked, not done. In Germany, for example, three quarters of the population prefer to retire before the age of 60, and a new immigration law is not even called by that name. Figure 3a also shows what would have happened if the birth rate would not have dropped, or if the characteristic age $`X`$ would not increase. Then, the increase in the fraction of old people would be much weaker. More realistically, Fig.3b shows the changes made by an immigration of people aged 6 to 40, amounting to half a percent of the total population per year and roughly stabilizing the population.

In summary, the possible breaks in the time development of demographic parameters were included in a simple approximation which could be incorporated into more detailed traditional models of human population dynamics. Simulations are in preparation  to check whether these phenomenological assumptions can be reproduced by more microscopic models .

I thank S. Moss de Oliveira for discussions.

 Wachter, K.W., Finch, C.E., 1997. Between Zeus and the Salmon. The Biodemography of Longevity, National Academy Press, Washington DC; see also La Recherche 322 (various authors), July/August 1999; Nature 408, No. 680 (various authors), November 9, 2000.

 Bomsdorf, E., 1993. Generationensterbetafeln für die Geburtsjahrgänge 1923-1993: Modellrechnungen für die Bundesrepublik Deutschland, Verlag Josef Eul, Köln.

 Tuljapurkar, S., Li, N., Boe, C., 2000. A universal pattern of mortality decline in the G7 countries. Nature 405, 789-792.

 Olshansky, S.J., Carnes B.A., Désesquelles A., 2001. Demography - Prospects for human longevity Science 291, 1491-1492.

 Wilmoth, J.R., Deegan, L.J., Lundström, H., Horiuchi, S., 2000. Increase of maximum life-span in Sweden, 1861-1999. Science 289, 2366-2368; Wilmoth, J.R., Horiuchi, S., 1999. Rectangularization revisited: Variability of age at death within human populations. Demography 36, 475-495; Wilmoth, J.R., Lundström, H., 1996. Extreme longevity in five countries - Presentation of trends with special attention to issues of data quality. Eur. J. Population 12, 63-93.

 Azbel, M. Ya., 1999. Phenomenological theory of mortality evolution: Its singularities, universality, and superuniversality. PNAS USA 96, 3303-3307.

 Azbel, M.Ya., 1999. Empirical laws of survival and evolution: Their universality and implications. PNAS USA 96, 15368-15373; 2001. Phenomenological theory of survival. Physica A 297, 235-241.

 Strehler, B. L., Mildvan, A.S. 1960. General theory of mortality and aging, Science 132, 14-21; L.A. Gavrilov, L.A., Gavrilova, N.S. 1991. The Biology of Life Span, Harwood Academic Publisher, Chur and 2001. The reliability theory of aging and longevity, J. Theor. Biology 213, 527-545; Azbel, M. Ya., 1996. Unitary mortality law and species-specific age. Proc. Roy. Soc. B 263, 1449-1454.

 Thatcher, A.R., 1999: The long-term pattern of adult mortality and the highest attained age. J. Roy. Statist. Soc. A 162, 5-30, and priv.comm.

 Thatcher, A.R., Kannisto, V., Vaupel, J.W., 1998. The Force of Mortality at Ages 80 to 120: Odense University Press, Odense. Robine, J.-M., Vaupel, J.W., 2001. Supercentenarians: slower ageing individuals or senile elderly?, Exp. Gerontology 36, 915-930 (2001)

 Finch C.E., Pike, M.C., 1996. Maximum life span predictions from the Gompertz mortality model. J. Gerontology (Biol. Sci.) 51 A, B 183-194.

 Yashin, A.I., Begun, A.S., Boiko, S.I., Ukraintseva, S.V., Oeppen, J., 2001. The new trends in survival improvement require a revision of traditional gerontological concepts., Exp. Gerontology 37, 157-167.

 Cebrat, S., 2002, private communication

 de Oliveira, P.M.C., Moss de Oliveira, S., Stauffer, D., Cebrat, S. 1999. Penna ageing model and improvement of medical care in 20th century. Physica A 273, 145-149; Niewczas, E., Cebrat, S., Stauffer, D. 2000. The Influence of the Medical Care on the Human Life Expectancy in 20th Century and the Penna Ageing Model. Theory Biosci. 119, 122-131; Moss de Oliveira, M., de Oliveira, P.M.C., Stauffer, D., 1999. Evolution, Money, War and Computers, Teubner, Stuttgart and Leipzig; Stauffer D., Radomski, J.R., 2001. Social effects in simple computer model of aging, Exp. Gerontology 37, 175-180.


format: text

text:
# Universality in escape from a modulated potential well

## Abstract

We show that the rate of activated escape $`W`$ from a periodically modulated potential displays scaling behavior versus modulation amplitude $`A`$. For adiabatic modulation of an optically trapped Brownian particle, measurements yield $`\mathrm{ln}W(A_\mathrm{c}A)^\mu `$ with $`\mu =1.5`$. The theory gives $`\mu =3/2`$ in the adiabatic limit and predicts a crossover to $`\mu =2`$ scaling as $`A`$ approaches the bifurcation point where the metastable state disappears.

Activated processes, such as escape from a metastable state, are exponentially sensitive to time-dependent external perturbations. Escape in driven systems has been investigated in studies of Josephson junctions Larkin ; Devoret  and infrared photoemission photoemission . It has been also explored recently in the context of thermally activated magnetization reversal in nanomagnets driven by a time-dependent magnetic field Wernsdorfer97  or a spin-polarized current Ralph02 . The strong dependence of the rates of activated processes on the wave form of the driving field enables selective control of the rates, which is important for many applications, for example control of rate and direction of diffusion in spatially periodic systems ratchets .

In this paper we investigate the dependence of the rate of activated escape $`W`$ on the amplitude $`A`$ of a periodic driving field. We show that this dependence displays universal behavior for large $`A`$ and investigate the scaling of $`\mathrm{ln}W`$ with $`A`$. Experimental results are obtained for a colloidal particle trapped in a modulated optically-generated potential. This system was used earlier McCann-99  for a quantitative test of the Kramers theory Kramers  of activated escape over a stationary potential barrier. We observe a power-law behavior of $`\mathrm{ln}W`$ as $`A`$ approaches a critical value. The critical exponent for adiabatic modulation is $`\mu =1.5`$, in agreement with theory. Adiabaticity is always violated close to the bifurcation point, and we predict the emergence of a different scaling, with $`\mu =2`$.

Field dependence of the escape rate is well understood for adiabatically slow modulation, where the field frequency $`\omega _F`$ is small compared to the reciprocal relaxation time $`t_\mathrm{r}^1`$. The driven system remains instantaneously in thermal equilibrium. The adiabatic escape rate $`W_{\mathrm{ad}}\mathrm{exp}[R_{\mathrm{ad}}/k_BT]`$ depends on the field primarily through the instantaneous activation energy $`R_{\mathrm{ad}}(t)`$. Even where the modulation of $`R_{\mathrm{ad}}`$ is small compared to the zero-field $`R_{\mathrm{ad}}`$, it may still substantially exceed $`k_BT`$, leading to strong modulation of $`W_{\mathrm{ad}}`$.

Finding the escape rate becomes more complicated for nonadiabatic modulation, where $`\omega _Ft_\mathrm{r}^1`$, because a driven system is far from thermal equilibrium. The problem has been addressed theoretically for different models of fluctuating systems Larkin ; Graham-84 ; Dykman-97ab ; SDG ; Hanggi-00 ; M&S-01 . The underlying idea is that activated escape results from an optimal fluctuation, which brings the system from a metastable to an appropriate saddle-type periodic state. We use it here to draw general conclusions about the scaling of the period-averaged escape rate $`\overline{W}`$ with the field amplitude $`A`$.

The field dependence of $`\mathrm{ln}\overline{W}`$ for $`\overline{W}t_\mathrm{r}^1,\omega _F`$ is sketched in Fig. 1. For small $`A`$, the field is a perturbation. Its major effect is to heat the system, with the change of the effective energy-dependent temperature quadratic in $`A`$ (the lowest order term after period averaging). As a result, $`\mathrm{ln}\left[\overline{W}/W_0\right]A^2`$, where $`W_0`$ is the escape rate in the absence of modulation Larkin ; Devoret . The range of $`A`$ is limited to $`\mathrm{ln}\left[\overline{W}/W_0\right]1`$.

For stronger modulating fields, the activation energy of escape $`R\mathrm{ln}\overline{W}`$ becomes linear in the amplitude $`A`$ Dykman-97ab . This behavior occurs because the optimal fluctuation leading to escape is a real-time instanton Langer , with duration $`t_\mathrm{r}`$. In the absence of modulation, escape has the same probability to happen at any time. Modulation lifts the time degeneracy and synchronizes escape events. As a result, even for zero-mean modulation $`\mathrm{ln}\overline{W}/W_0A`$. For an overdamped Brownian particle, an explicit solution was obtained SDG  throughout the heating and log-linear regions in Fig. 1, and the results have been confirmed by extensive simulations Chaos-01 .

As we show in this paper, there is another region where the dependence of $`\overline{W}`$ on $`A`$ is universal and unexpected: the vicinity of the saddle-node bifurcation point $`𝒜_\mathrm{c}`$ where the metastable state of forced vibrations merges with a saddle-type periodic state Guckenheimer . In this range one of the motions of the system becomes slow. The universality is related to the corresponding critical slowing down. However, it turns out that the system can display different types of critical behavior depending on the relationship between $`\omega _F`$ and the relaxation time in the absence of driving $`t_{\mathrm{r0}}`$.

Experiments on amplitude-dependent escape from a sinusoidally modulated potential were performed using an optically-trapped Brownian particle in water. A configurable optical potential was constructed by focusing two parallel laser beams through a single microscope objective lens. Each beam creates a stable three-dimensional trap as a result of electric field gradient forces exerted on a transparent dielectric spherical silica particle of diameter 0.6 $`\mu `$m. At the focal plane, the beams are typically displaced by 0.25 to 0.45 $`\mu `$m, creating a double-well potential. The depths of the potential wells and the height of the intervening barrier, which is typically 1 to 10 $`k_BT`$, are readily controlled by adjusting beam intensities and separations.

The two HeNe lasers (17 mW, 633 nm) that create the traps are stabilized by electro-optic modulators and imaged into a sample cell by a 100x objective. The beams are mutually incoherent and circularly polarized when they enter the microscope. A single trapped colloidal sphere is imaged onto a digital camera operating at 200 frames/s. The $`x`$ and $`y`$ coordinates of the sphere’s center are computed using a pattern matching routine that yields spatial resolution better than 10 nm, whereas the $`z`$ coordinate parallel to the light wave vector is extracted by analysis of the image as the particle position fluctuates about the focal plane. The overall analysis and storing of the sphere’s $`(x,y,z)`$ coordinates is accomplished in less than a frame duration so that no images need to be recorded.

The stability of the system is sufficient to compute the full 3-dimensional optical potential from long-time measurements of the particle probability density McCann-99 . The absolute transition rates calculated from the Kramers expression using measured curvatures at the stable and unstable points of the potential are in quantitative agreement with the experimental data over 4 orders of magnitude.

The electro-optic modulators that stabilize the optical traps also enable ac modulation with arbitrary waveforms derived from an electronic signal generator. When the two beams overlap at the separations required to form a potential barrier in the range 5 to 8 $`k_BT`$, each contributes to the position and shape of both potential minima and the potential in the intermediate range. As a consequence of this nonlocality, a small modulation of one beam causes a nearly antisymmetric change of both barrier heights.

The results on over-barrier transitions from one of the wells for the modulating signal $`F(t)=A\mathrm{cos}\omega _Ft`$ are plotted in Fig. 2. The amplitude $`A`$ in Fig. 2 is normalized to $`k_BT`$, where T = 298 K. The intrawell relaxation rate in the absence of modulation $`t_{\mathrm{r0}}^110^3`$ s<sup>-1</sup> McCann-99  is much larger than $`\omega _F`$, so the data fall into the adiabatic regime. In this case escape is most likely to occur when the instantaneous potential barrier $`\mathrm{\Delta }U(t)`$ is minimal. For large modulation amplitude and for temperatures used here, the maximal escape rate $`W_m`$ becomes comparable to the modulation frequency, and therefore it is $`W_m`$ rather than the period-averaged rate $`\overline{W}`$, that is the more relevant measure. In obtaining $`W_m`$ we used the distribution of the dwell time illustrated in the inset to Fig. 2. The dwell time was computed for the modulated double-well system as the time interval between the appearance in, and escape from, a well.

The solid line in Fig. 2 represents a least-squares fit to the data for large $`A`$ using the function $`\mathrm{ln}W_m/W_0=C_1C_2(A_\mathrm{c}A)^\mu `$, where $`\mu =1.5`$, and the critical modulation amplitude $`A_\mathrm{c}`$ and $`C_{1,2}`$ are fitting parameters uncertainties .

To interpret these results we will consider activated escape in the system for $`A`$ close to the critical amplitude $`𝒜_\mathrm{c}`$ where the metastable and unstable periodic states merge. In the adiabatic approximation, the stable and unstable states correspond to the potential well and the barrier top. Once per period, when the driving $`|F(t)|`$ is a maximum (for example when $`t=n\tau _F,n=0,\pm 1,\mathrm{}`$), the stable and unstable states are closest to each other and the barrier height $`\mathrm{\Delta }U(t)`$ is a minimum. We take $`F(0)`$ positive, $`F(0)=A`$. The well and barrier top merge for $`A=A_\mathrm{c}`$, which is the adiabatic bifurcation point.

For small $`\delta A=AA_\mathrm{c}`$ and $`|tn\tau _F|`$, the dynamics of the system near the metastable state is slow Guckenheimer . It is described by the “soft mode” with coordinate $`q(t)`$. For negative $`\delta A`$ the system has a stable and an unstable state, $`q_a`$ and $`q_b`$. We set $`q=0`$ at the saddle-node state where they merge in the limit of small $`\omega _F`$ for $`A=A_\mathrm{c},t=n\tau _F`$. Close to $`A_\mathrm{c}`$ and $`n\tau _F`$ the equation of motion for $`q`$, neglecting fluctuations, takes the form

$$t_{\mathrm{r0}}\frac{dq}{dt}=q^2+F(t)A_\mathrm{c}$$
(1)
($`t_{\mathrm{r0}}`$ is the relaxation time for $`F=0`$, $`\omega _Ft_{\mathrm{r0}}1`$). In the adiabatic approximation, $`F(t)`$ depends parametrically on $`t`$. Activated transitions occur over the instantaneous barrier $`\mathrm{\Delta }U(t)=(4/3)[A_\mathrm{c}F(t)]^{3/2}`$.

The relaxation time of the system $`t_\mathrm{r}=t_{\mathrm{r0}}[A_\mathrm{c}F(t)]^{1/2}`$ becomes large close to the bifurcation point. Therefore even a slowly varying field will ultimately become too fast for the system to follow adiabatically. To allow for nonadiabaticity we expand $`F(t)`$ in the region of small $`\omega _F|tn\tau _F|`$. The equation of motion (1) in scaled variables $`Q=q/\mathrm{\Omega },\tau =\mathrm{\Omega }(tn\tau _F)/t_{\mathrm{r0}}`$ is then

$`\dot{Q}=K(Q,\tau ),K=Q^2+\zeta \tau ^2,`$ (2)
with $`\zeta =\mathrm{\Omega }^2(AA_\mathrm{c})`$. In Eq. (2) $`\dot{Q}dQ/d\tau `$, and $`\mathrm{\Omega }=(\frac{1}{2}t_{\mathrm{r0}}^2d^2F/dt^2)^{1/4}`$ (the derivative is evaluated at $`t=0`$ and $`A=A_\mathrm{c}`$). For sinusoidal modulation $`\mathrm{\Omega }=(A_\mathrm{c}t_{\mathrm{r0}}^2\omega _F^2/2)^{1/4}`$.

The dynamics (2) is determined by one dimensionless parameter $`\zeta `$. For $`\zeta 1`$, the adiabatic approximation applies, and the stable and unstable states are given by $`Q_{a,b}(\tau )=(\zeta +\tau ^2)^{1/2}`$. On the other hand, for $`|\zeta |1`$ the time dependence of the states becomes distorted and asymmetric, cf. inset in Fig. 3. For $`\zeta =1`$, the states given by Eq. (2) merge and $`Q_a(\tau )=Q_b(\tau )=\tau `$. The corresponding nonadiabatic bifurcational value of the amplitude is $`A_\mathrm{b}=A_\mathrm{c}+\mathrm{\Omega }^2`$.

Escape from the attractor can be quite generally described by adding a random Gaussian force to Eq. (1). Because the system dynamics is slow, this force is effectively $`\delta `$-correlated (white noise), with intensity $`Dk_BT`$ for thermal noise. The corresponding force $`\xi (\tau )`$ in Eq. (2) for $`\dot{Q}`$ has a correlator $`\xi (\tau )\xi (\tau ^{})=2\stackrel{~}{D}\delta (\tau \tau ^{}),\stackrel{~}{D}=D/\mathrm{\Omega }^3t_{\mathrm{r0}}.`$

The noise-driven system is most likely to escape during the time when $`F(t)`$ is close to its maximum. The total probability to escape over this time is $`W_t=\overline{W}\tau _F`$, where $`\overline{W}`$ is the period-averaged escape rate. We will calculate $`W_t`$ assuming that the noise intensity $`D`$ is the smallest parameter of the theory, in which case $`W_t\mathrm{exp}(R/D)`$. The calculation can be done by generalizing the instanton technique Chaos-01  to the Langevin equation $`\dot{Q}=K(Q,\tau )+\xi (\tau )`$. The activation energy of escape $`R=(D/\stackrel{~}{D})S`$ is given by the solution of the variational problem

$$S=\frac{1}{4}\mathrm{min}_{\mathrm{}}^{\mathrm{}}𝑑\tau (\dot{Q}K)^2.$$
(3)
The boundary conditions are $`Q(\tau )Q_a(\tau )`$ for $`\tau \mathrm{}`$ and $`Q(\tau )Q_b(\tau )`$ for $`\tau \mathrm{}`$. They correspond to the picture in which, prior to escape, the system is performing small fluctuations about the metastable state $`Q_a`$. Escape occurs as a result of a large fluctuation that brings the system to the saddle state $`Q_b`$. The most probable among such fluctuations is the one in which the system moves along a trajectory \[the most probable escape path (MPEP)\] that minimizes the functional (3). From Eq. (3), the activation energy $`R`$ is a function of one parameter $`\zeta `$, shown in Fig. 3.

An explicit expression for $`R`$ can be obtained in the adiabatic limit $`\zeta 1`$. Here, one can ignore the term $`\tau ^2`$ in $`K`$ (2) for typical $`|\tau |\zeta ^{1/2}`$. Then the MPEP is an instanton $`Q_{\mathrm{opt}}(\tau ,\tau _0)=|\zeta |^{1/2}\mathrm{tanh}\left[|\zeta |^{1/2}(\tau \tau _0)\right]`$ centered at an arbitrary $`\tau _0`$. The adiabatic activation energy is $`RS=(4/3)(\zeta )^{3/2}(A_\mathrm{c}A)^{3/2}`$ 3/2 . This agrees with the experimental data shown in Fig. 2.

The time-dependent term in $`K`$ lifts translational invariance of the instanton $`Q_{\mathrm{opt}}(\tau ,\tau _0)`$ and synchronizes escape events. The action $`S`$ is minimal for $`\tau _0=0`$. To lowest order in $`1/\zeta `$

$$S(4/3)|\zeta |^{3/2}+(\pi ^2/6)|\zeta |^{1/2}.$$
(4)
From (4), the nonadiabatic correction to the activation energy diverges as $`|\delta A|^{1/2}`$ as $`A`$ approaches $`A_\mathrm{c}`$.

The activation energy may also be found explicitly close to the bifurcation point $`\zeta =1`$. In this case escape is most likely to occur while the coexisting attracting and repelling trajectories $`Q_a(\tau )`$ and $`Q_b(\tau )`$ stay close to each other. Because of the special structure of the states $`Q_{a,b}(\tau )`$, the variational problem (3) for the MPEP can be linearized and solved, giving

$$S=(\pi /8)^{1/2}(1\zeta )^2,\mathrm{\hspace{0.33em}1}\zeta 1.$$
(5)
Eq. (5) shows that, near the nonadiabatic bifurcation point $`A_\mathrm{b}`$, there emerges another scaling region, where the activation energy $`R(A_\mathrm{b}A)^2`$. The nonadiabatic exponent is thus $`\mu =2`$, in contrast to $`\mu =3/2`$ in the adiabatic limit. It describes $`R(A)`$ for the modulation amplitude lying between $`A_\mathrm{c}`$ and $`A_\mathrm{b}`$, as seen from Fig. 3. The width of this interval is $`\omega _F`$ provided $`\omega _Ft_{\mathrm{r0}}1`$.

The above theory does not apply in the exponentially narrow range $`1\zeta \mathrm{exp}[C/\omega _Ft_{\mathrm{r0}}]`$ ($`C1`$). This is because, as $`A`$ approaches $`A_\mathrm{b}`$, the relaxation time of the system $`t_\mathrm{r}`$ becomes logarithmically long. Then the expansion of the driving force in $`t/\tau _F`$ that leads to the equation of motion (2) becomes inapplicable. For $`\omega _Ft_\mathrm{r}1`$ the duration of the MPEP greatly exceeds the period $`\tau _F`$. In this case the activation energy again displays the $`\mu =3/2`$ power-law dependence on the distance to the true bifurcation value of the amplitude $`𝒜_\mathrm{c}`$. The $`(𝒜_\mathrm{c}A)^{3/2}`$-law also applies asymptotically if the driving is nonadiabatic for small $`A`$.

The experimental data discussed above refer to very small $`\omega _Ft_{\mathrm{r0}}<10^3`$. Therefore the nonadiabatic region of modulating amplitudes is narrow and the activation energy there is less than $`k_BT`$. Higher modulation frequencies will be required to detect the crossover from the observed $`\mu =3/2`$ scaling to the nonadiabatic $`\mu =2`$ scaling. It appears that the $`\mu =3/2`$ scaling function describes the data in a broad range of $`A`$. It thus provides a good interpolation of the activation energy from the linear in $`A`$ to the critical region.

In conclusion, we have identified several regions where the activation energy $`R`$ of escape from a metastable potential displays scaling behavior as a function of the amplitude $`A`$ of an externally applied periodic modulation. We have demonstrated with a system whose potential was directly measured that, near the amplitude $`A_\mathrm{c}`$ where the local minimum and maximum of the potential contact one another, $`R(A_\mathrm{c}A)^{3/2}`$. We have also shown that, because of nonadiabaticity near the bifurcation point, there necessarily emerges a region where the scaling exponent changes from $`\mu =3/2`$ to $`\mu =2`$. We expect that these scalings can be observed in other systems, such as modulated Josephson junctions Devoret  and nanomagnets Wernsdorfer97 ; Ralph02 ; Koch00 . The ideas discussed here can be used to gain additional insight into the physics of these systems and more generally, activated effects in nonequilibrium systems.

This research was supported by the NSF DMR-9971537 and NSF PHY-0071059.


text:
# CERN-TH/2002-077 hep-ph/0204051 April, 2002 The Problem of Large Leptonic Mixing

## Acknowledgements

I am grateful to G.C. Branco for a thorough reading of the manuscript and suggestions. This work received partial support from the Portuguese Ministry of Science - Fundação para a Ciência e Tecnologia and the CERN Theory Division which I thank for the kind hospitality extended.


text:
# REALISTIC INTERATOMIC POTENTIAL FOR MD SIMULATIONS

## 1. Introduction

Investigation of strongly anharmonic nonlinear atomic systems by molecular dynamics (MD) method at high temperatures, pressures, or study of systems effected by large amplitude excitations requires high accuracy of interatomic potential (IP). Series expansion of the IP in the displacements of atoms from equilibrium positions is widely used both in phonon theory and in MD simulation . Usually, fourth-order anharmonisms or lower-order ones can be taken into account because of complexity of expansion coefficients calculation. As an alternative, realistic potential method is used , in which exact equations of motion of atoms are solved using IP of concrete substance without series expansion. Owing to that, all-order anharmonisms are taken into account automatically. This advantage of realistic potential method is especially useful in the MD simulation of soliton solutions where atoms approach each other closely. Realistic IP should have simplest form to reduce calculation expenses as well as it must describe precisely the properties of the substance under extreme conditions. The aim of this paper is to obtain such IP.

Conventional way of realistic IP determination is empirical fitting to the properties of gas or a crystal near the equilibrium point . However, such potentials become unreliable at small interatomic distances like to that arising in soliton waves. The properties of highly compressed matter (e.g., for neon up to 1 Mbar ) could give an information for obtaining all-distance reliable IP. However, the set of properties, which can be measured accurately at megabar pressures is restricted strongly. Practically, the equation of state and bulk modulus only may be included in this set . For $`C_{ik}`$ modules the precision worsens drastically even at kilobar pressures . There is insufficiency of empirical information for fitting all the parameters of IP, and ab-initio calculation is required.

Realistic IP via interatomic distance is obtained in present work for the crystal and dimer of neon. Repulsive part of the potential is calculated ab-initio in Hartree-Fock approximation using the basis of localized atomic orbitals orthogonalized exactly on different lattice sites. Attractive part is chosen to have standard Van-der-Vaals form of $`Cr^6`$ with single empirical parameter $`C`$. Used approximations and calculation details are described in the section 2. In the section 3 repulsive part of IP is interpolated by exponential function of interatomic distance (Exp-6 potential) and the IP parameter are determined. Experimental verification of the IP obtained is performed in the section 4 using the data concerning equation of state  and elastic moduli  of compressed neon. The IP calculated is found to be in a good consistence with the experiment in whole range of pressure.

## 2. Ab-initio calculation of repulsion potential

In MD simulations by realistic potential method the problem is divided into two stages. The former is quantum-mechanical calculation of the IP at electron level, with interatomic distance considering as a parameter. The latter is solving equations of motion of atoms using the IP obtained. This division is correspondent to adiabatic approximation when atoms and electrons motion is described separately .

Since pair collisions of atoms have maximal probability, we concentrate the attention on the dimer of neon, and define the IP as a cohesive energy of the dimer. Three-atom forces can be taken into account as a correction to the two-atom ones using incremental expansion . The estimation of  shows three-atom force contribution to be small.

In Hartree-Fock approximation short-range repulsive part of IP is expressed through one-electron density matrix. We don’t use hard core approximation. Rearrangement of all electron shells is allowed as interatomic distance is altered.

Localized basis of atomic orbitalls orthogonalized exactly (by Lovdin procedure ) on different lattice sites is used. In this basis one-electron density matrix has the form

$`\rho (\text{r}^{}|\text{r};\{\text{l}\})=2{\displaystyle \underset{\text{l}s}{}}\{\phi _s(\text{r}^{}\text{l})\phi _s^{}(\text{r}\text{l}){\displaystyle \underset{\text{l}^{}s^{}}{}}\phi _s^{}(\text{r}^{}\text{l}^{})P_{s^{}s}^{\text{l}^{}\text{l}}\phi _s^{}(\text{r}\text{l})\},`$
(2.1) $`\text{P}=\text{I}(\text{I}+\text{S})^1,`$
where $`\phi _s(\text{r}\text{l})=|\text{l}s>`$ is wave function of electron in isolated atom (atomic orbital), l and $`\text{l}^{}`$ are radius-vectors of lattice sites, $`s`$ numerates occupied states of the atom, P is orthogonalizing matrix, I is unit matrix, S is overlap integral matrix with the elements

$`S_{s^{}s}^{\text{l}^{}\text{l}}=<\text{l}^{}s^{}|\text{l}s>;\text{l}\text{l}^{},`$
(2.2) $`S_{s^{}s}^{\text{l}^{}\text{l}}=0;\text{l}=\text{l}^{}.`$
We expand repulsive part of IP in the terms of small parameter such as the largest overlap integral $`S`$. Usually, $`S<<1`$ in uncompressed crystal, and overlap integrals grow exponentially as interatomic distance is decreased. The IP is expressed through the products of elements of orthogonalizing matrix P and two-center Slater-Koster integrals. These integrals are atomic obital matrix elements of crystal hamiltoinial operators. The order in $`S`$ for two-center integrals is estimated using the theorem about average value. The elements of matrix $`\text{P}=\text{I}(\text{I}+\text{S})^1`$ are expanded in powers of overlap integrals matrix S

$`P_{ss^{}}^{\text{l}\text{l}^{}}=S_{ss^{}}^{\text{l}\text{l}^{}}+O(\text{S}^2),`$
(2.3) $`P_{ss^{}}^{\text{l}\text{l}}=(\text{S}^2)_{ss^{}}^{\text{l}\text{l}}+O(\text{S}^3).`$
The elements of matrix P contain high-order terms along with the main ones proportional to $`S`$ and $`S^2`$.

Using the estimations described above, we expand the repulsive part of IP in powers of $`S`$

(2.4) $`V_{sr}=E^{(0)}+W_2+W_4+W_6.`$
Here $`E^{(0)}`$ is the energy of interatomis interaction if orthogonalizing of neighbor atoms orbitals is neglected, $`W_2,W_4,W_6`$ are orthogonalizing corrections. Series expansion in $`S`$ begins for them from the second, the third, and the sixth powers respectively. Due to the presence of matrix P, orthogonalizing corrections contain high-order terms in $`S`$ along with the main ones.

In the equation 2.4

(2.5) $`E^{(0)}={\displaystyle \underset{\text{l}s}{}}{\displaystyle \underset{\text{m},\text{m}\text{l}}{}}\text{l}s|V_{en}^\text{m}+V_a^\text{m}+V_{ex}^\text{m}|\text{l}s+U_{nn}.`$
The first term in equation 2.5 consists of two-center integrals. They are atomic orbital matrix elements of electron-ion interaction potential $`V_{en}^\text{m}`$, of neutral isolated atom potential $`V_a^\text{m}`$, of electron-electron exchange interaction potential $`V_{ex}^\text{m}`$ respectively. The second term is the energy of nucleus-nucleus interaction. Electron-ion interaction potential has the form

(2.6) $`V_{en}^\text{m}=V_{en}(\text{r}\text{m})=Ze^2/|\text{r}\text{m}|.`$
Neutral isolated atom potential is

(2.7) $`V_a^\text{m}=V_a(\text{r}\text{m})=V_{en}(\text{r}\text{m})+2{\displaystyle \underset{t}{}}<\text{m}t|v_c|\text{m}t>,`$
where

$`<\text{m}t|v_c|\text{m}t>={\displaystyle \phi _t^{}(\text{r}^{}\text{m})v_c(\text{r}\text{r}^{})\phi _t(\text{r}^{}\text{m})𝑑\text{r}^{}},`$
$`v_c(\text{r}\text{r}^{})=e^2/|\text{r}\text{r}^{}|.`$
Action of electron-electron exchange interaction potential on wave function is defined as

(2.8) $`<\text{l}s|V_{ex}^\text{m}|\text{l}s>={\displaystyle \underset{t}{}}<\text{l}s,\text{m}t|v_c|\text{l}s,\text{m}t>.`$
In the equation 2.4 orthogonalizing corrections, $`W_2,W_4,W_6`$, are of the form

$`W_2=2{\displaystyle \underset{\text{l}s}{}}{\displaystyle \underset{\text{l}^{}s^{},\text{l}\text{l}^{}}{}}P_{ss^{}}^{\text{l}\text{l}^{}}<\text{l}^{}s^{}|V_a^\text{l}^{}+V_{ex}^\text{l}^{}|\text{l}s>`$
(2.9) $`{\displaystyle \underset{\text{lm}ss^{}tt^{},\text{l}\text{m}}{}}P_{ss^{}}^{\text{ml}}P_{tt^{}}^{\text{lm}}<\text{l}s^{},\text{m}t^{}|v_c|\text{m}s,\text{l}t>;`$
$`W_4={\displaystyle \underset{\text{l}ss^{}tt^{}}{}}P_{ss^{}}^{\text{l}\text{l}}P_{tt^{}}^{\text{l}\text{l}}\{2<\text{l}s^{},\text{l}t^{}|v_c|\text{l}t,\text{l}s><\text{l}s^{},\text{l}t^{}|v_c|\text{l}s,\text{l}t>\}`$
$`2{\displaystyle \underset{\text{l}ss^{}}{}}P_{ss^{}}^{\text{l}\text{l}}<\text{l}s^{}|{\displaystyle \underset{\text{m}\text{l}}{}}(V_a^\text{m}+V_{ex}^\text{m})|\text{l}s>+`$
$`+2{\displaystyle \underset{\text{l}\text{m}ss^{}tt^{},\text{l}\text{m}}{}}\{P_{ss^{}}^{\text{l}\text{l}}P_{tt^{}}^{\text{m}\text{m}}<\text{l}s^{},\text{m}t^{}|v_c|\text{m}t,\text{l}s>+P_{ss^{}}^{\text{m}\text{l}}P_{tt^{}}^{\text{l}\text{m}}<\text{l}s^{},\text{m}t^{}|v_c|\text{l}t,\text{m}s>\}+`$
$`+2{\displaystyle \underset{\text{l}\text{m}ss^{}tt^{},\text{l}\text{m}}{}}P_{ss^{}}^{\text{m}\text{l}}P_{tt^{}}^{\text{m}\text{l}}\{2<\text{l}s^{},\text{l}t^{}|v_c|\text{m}t,\text{m}s><\text{l}s^{},\text{l}t^{}|v_c|\text{m}s,\text{m}t>\}+`$
(2.10) $`+4{\displaystyle \underset{\text{l}\text{m}ss^{}tt^{},\text{l}\text{m}}{}}P_{ss^{}}^{\text{m}\text{l}}P_{tt^{}}^{\text{l}\text{l}}\{2<\text{l}s^{},\text{l}t^{}|v_c|\text{l}t,\text{m}s><\text{l}s^{},\text{l}t^{}|v_c|\text{m}s,\text{l}t>\};`$
(2.11) $`W_6={\displaystyle \underset{\text{l}\text{m}ss^{}tt^{},\text{l}\text{m}}{}}P_{ss^{}}^{\text{l}\text{l}}P_{tt^{}}^{\text{m}\text{m}}<\text{l}s^{},\text{m}t^{}|v_c|\text{l}s,\text{m}t>.`$
Since the orthogonalizing corrections grow exponentially as the interatomic distance is decreased it is impossible to say what correction may be neglected. It should be checked for each substance under consideration.

Using the method described, we calculate repulsive part of IP, $`V_{sr}`$ (equation 2.4), for neon dimer as a function of interatomic distance $`d`$. Atomic orbitals from Clementi-Roetti set  are used as a basis. Hartree system of atomic units $`\mathrm{}=e=m_e=1`$ is applied. The calculation shows the terms $`E^{(0)}`$ and $`W_2`$ in equation 2.4 to have the same order of magnitude and opposite signs. These terms are found to give major contributions to the IP. The $`W_4`$ correction consists of 0,02 per cent of the IP at equilibrium interatomic distance $`d_0`$. Further, the $`W_4`$ does not exceed of 1 per cent of the IP up to $`d0.75d_0`$. Finally, at small $`d`$, like to that arising in soliton waves (for $`d`$ above 0,6-0.75 $`d_0`$), the $`W_4`$ becomes about 2-4 per cent of the IP. The contribution of $`W_6`$ to the IP is small negligibly (0.002 per cent) in whole range of $`d`$ under consideration.

## 3. Determination of interatomic potential parameters

We interpolate calculated points $`V_{sr}(d)`$ by exponential function of interatonic distance using least square method by the formula

$`V_{sr}(d)=A_0\mathrm{exp}(\alpha (x1));`$
(3.1) $`x=d/z_0`$
with two unknown parameters $`A_0`$ and $`\alpha `$. Experimental equilibrium interatomic distance for neon dimer $`z_0=5,8411`$ a.u.  is used as the third parameter of the IP. The parameters are found to be $`A_0=(1,1384\pm 0,0002)10^4`$ a.u., $`\alpha =13,6407\pm 0,0037`$. Interpolation error is 4-1 per cent of $`V_{sr}`$ when the $`d`$ is altered from equilibrium one to 0.6$`z_0`$.

Adding the attractive part, we express the IP in standard Exp-6 form

$`V(d)=A_0\mathrm{exp}(\alpha (x1))Cd^6;`$
(3.2) $`x=d/z_0.`$
A single unknown parameter $`C`$ remains in attractive part of IP. We propose to fit the $`C`$ to experimental equilibrium interatomic distance. Using of equilibrium data is considered to be reliable at all interatomic distances since the attraction is essential near the equilibrium only while ab-initio calculated repulsive part dominates at small $`d`$.

For MD simulation of lattice dynamics, it is possible to fit the $`C`$ to experimental data for dimer at $`T`$=0 K because the temperature effects will be taken into account explicitly, at the stage of equations of motion solving. In this case, for neon $`C`$=10,7293 (experimental equilibrium interatomis distance in the dimer is $`z_0`$=5,8411 a.u. ). Calculated cohesive energy of dimer is $`E_{coh}=1,449710^4`$ a.u., experimental one is $`E_{coh}=1,33810^4`$ a.u. . The discrepancy is 7 per cent of experimental value.

For calculating static properties of a crystal at finite temperature, e.g., equation of state, elastic modules, it is better to fit the $`C`$ to experimental data for a crystal at the same temperature. Such determination allows one to take into account implicitly three-atom forces, temperature effects, zero-point oscillations, and other effects omitted at the stage of IP calculating. In this case, for neon $`C`$=7,4030 (experimental equilibrium interatomis distance in the crystal is $`d_0`$=5,9647 a.u. at $`T=4,25`$ K ). Calculated cohesive energy of uncompressed crystal is $`E_{coh}=6,762010^4`$ a.u. per atom, experimental one is $`E_{coh}=(7,35\pm 0,03)10^4`$ a.u. . The discrepancy is 7.6 per cent of experimental value.

## 4. Results and discussion

Interatomic potential of neon is given in the figure 1 as a function of interatomic distance $`d`$. The IP calculated by equation 3 for dimer is plotted by solid curve. Van-der-Vaals constant ($`C`$=10,7293) is fitted to experimental equilibrium interatomic distance in dimer .

”Experimental” IP obtained in  is denoted by solid circles. This IP had been determined by interpolating experimental data $`p(V)`$ (measured at 300 K) by the formula Exp-6. The interpolation had been performed in theoretical model taking thermal pressure and zero-point oscillations into account explicitly, excluding them from the definition of IP. It allows us to compare the 300 K data of  with our zero-temperature result. Three-atom forces didn’t include explicitly in the model of . However, in , the effect of these forces is taken into account implicitly through fitting the IP to experimental data for a crystal. In our calculation three-atom forces are omitted because of fitting to dimer data. The agreement of calculated IP and experimental one indicates that three-atom forces in neon are small at the pressures up to 1Mbar.

Two remaining curves in the figure 1 are interatomic potentials of neon obtained by fitting to experimental data using Lennard-Jones potential (6-12 formula)

$`V(x)=\epsilon (2/x^6+1/x^{12})`$
$`x=d/z_0,`$
where $`\epsilon `$ and $`z_0`$ are fitting parameters. Dashed curve is the IP obtained using corresponding-states law fitted to vapor-pressure ratio of isotopic liquid . Dashed-dotted curve is the IP fitted to experimental lattice constant and cohesive energy of crystal neon at $`p=0`$, $`T=0`$ K . Fitting to equilibrium crystal properties leads to bad describing the IP for compressed crystal. Fitting to compressed gas properties gives the values of the IP close to experimental ones at moderated pressures.

Using the IP obtained (eq. 3) we calculate the equation of state $`p(V)`$ for solid neon. Calculated pressure p against fractional volume is given in the figure 2 as solid curve. Van-der-Vaals constant is fitted to experimental equilibrium interatomic distance $`d_0`$=5,9647 a.u. measured for crystal neon at $`T`$=4.25 K, $`p`$=0 .

Experimental points $`p(V)`$ from  ($`T`$=300 K) and  ($`T`$=4.2 K) are also given in the figure 2. At the pressures below 20 kbar theoretical curve is in a good agreement with the experimental points of . At moderated pressures theoretical curve deviates from experimental points of  by 4 per cent. This deviation caused, mainly, by neglecting of thermal pressure in our calculation. Figure 2 shows temperature sensitivity of the equation of state to be small.

We calculate bulk modulus of solid neon by means of the IP obtained. Van-der-Vaals constant is fitted to experimental equilibrium interatomic distance in the crystal . Calculated bulk modulus $`B`$ via the pressure $`p`$ is given in the figure 3 as solid curve. Experimental points obtained in  at $`T`$=4.2 K are plotted as solid symbols. Bulk modulus is seen to be more sensitive to the approximations used. Growing when the $`p`$ is enhanced, the difference between calculated $`B`$ and measured one becomes about 7 per cent of experimental $`B`$ at $`p`$=20 kbar. Incorrect taking three-atom forces into account at moderated pressures is seems to contribute mainly in this discrepancy. In our calculation three-atom forces (and zero-point oscillations too) are taken into account implicitly, by fitting the IP to experimental data for uncompressed crystal. Thus, calculated $`B`$ agrees with experimental one at small pressures only (to 8 kbar). One can’t determine correctly the dynamics of alteration of three-atom forces with enhancing of pressure. It is the cause of growing the deviation of calculated $`B`$ from measured one.

We calculate elastic modules $`C_{ik}`$ using the IP obtained with Van-der-Waals constant fitted to crystal experimental data . Calculated modules and experimental ones are given in table 1 for uncompressed solid neon at low temperatures. Isothermic modules had been obtained in static measurements . Adiabatic modules had been measured in ultrasonic and neutron scattering experiments . However, the difference between isotermic modules and adiabatic ones is negligible at the temperatures under consideration (see, e.g., ).

The $`C_{ik}`$ modules are seen to be more sensitive to the measurement method and calculation approximations. The difference between theoretical and experimental values of $`C_{11}`$ and $`C_{44}`$ is about 10 per cent of experimental values for most accurate experiment . The agreement is better for $`C_{12}`$ modulus (the discrepancy is about 2 per cent ). The deviation from Cauchy relation $`\delta =(C_{44}C_{12})/C_{12}`$ is also given in table 1. Cauchy violation is the measure of deviation of the IP from spherical symmetry. The $`\delta =0,11\pm 0,03`$ in , while it falls into experimental error bar in other experiments listed in the table 1. Cauchy relation takes place for our calculation results because spherical symmetry form of the IP is supposed in theoretical model. Small value of experimental $`d`$ indicates that spherical symmetry approximation for IP is valid for uncompressed neon at least. For another rare gas crystal, krypton, experiment  shows Cauchy relation to satisfy well under pressure up to 80 kbar. Moreover, for MgO the Cauchy violation is measured to drop with enhancing pressure up to 200 kbar .

Unlike to $`C_{ik}`$ modules, bulk modulus $`B`$ is less sensitive to measurement method and calculation approximations. The discrepancy of theoretical result and experimental one doesn’t exceed of 4 per cent and falls into experimental error frames.

## 5. Conclusion

Coefficients of realistic IP of simple form Exp-6 are obtained for neon by ab-initio calculation of repulsive part in Hartree-Fock approximation in the basis of atomic orbitals orthogonalized exactly on different lattice sites. Attractive part is determined empirically using single fitting parameter, Van-der-Vaals constant $`C`$. For fitting the $`C`$ it is enough to know experimental equilibrium interatomic distance in crystal (or dimer), i.e. high pressure experimental data is not required. The IP calculated is suitable for molecular dynamic simulations of high temperature and high pressure properties and phenomena in crystals and liquids due to simplicity of the form and precise describing experimental data in wide range of interatomic distances.