text:
# Power-law persistence and trends in the atmosphere: A detailed study of long temperature records

## I Introduction

The persistence of weather states on short terms is a well-known phenomenon: A warm day is more likely to be followed by a warm day than by a cold day and vice versa. The trivial forecast, that the weather of tomorrow is the same as the weather of today, was in previous times often used as a “minimum skill” forecast for assessing the usefulness of short-term weather forecasts. The typical time scale for weather changes is about 1 week, a time period that corresponds to the average duration of so-called “general weather regimes” or “Grosswetterlagen”, so this type of short-term persistence usually stops after about 1 week. On larger scales, other types of persistence occur. One of them is related to circulation patterns associated with blocking . A blocking situation occurs when a very stable high pressure system is established over a particular region and remains in place for several weeks. As a result the weather in the region of the high remains fairly persistent throughout this period. It has been argued recently  that this short-term persistence regime may be linked to solar flare intermittency. Furthermore, transient low pressure systems are deflected around the blocking high so that the region downstream of the high experiences a larger than usual number of storms. On even longer terms, a source for weather persistence might be slowly varying external (boundary) forcing such as sea surface temperatures and anomaly patterns. On the scale of months to seasons, one of the most pronounced phenomena is the El Nino southern oscillation event which occurs every 3–5 years and which strongly affects the weather over the tropical Pacific as well as over North America .

The question is, how the persistence that might be generated by very different mechanisms on different time scales decays with time $`s`$. The answer to this question is not easy. Correlations, and in particular long-term correlations, can be masked by trends that are generated, e.g., by the well-known urban warming. Even uncorrelated data in the presence of long-term trends may look like correlated ones, and, on the other hand, long-term correlated data may look like uncorrelated data influenced by a trend.

Therefore, in order to distinguish between trends and correlations one needs methods that can systematically eliminate trends. Those methods are available now: both wavelet techniques (WT) (see, e.g., Refs. ) and detrended fluctuation analysis (DFA) (see, e.g., Refs. ) can systematically eliminate trends in the data and thus reveal intrinsic dynamical properties such as distributions, scaling and long-range correlations very often masked by nonstationarities.

In a previous study , we have used DFA and WT to study temperature correlations in different climatic zones on the globe. The analysis focused on 14 continental stations, several of them were located along coastlines. The results indicated that the temperature variations are long-range power-law correlated above some crossover time that is of the order of 10 days. Above the crossover time, the persistence, characterized by the autocorrelation $`C(s)`$ of temperature variations separated by $`s`$ days, decayed as

$$C(s)s^\gamma ,$$
$`(1)`$
where, most interestingly, the exponent $`\gamma `$ had roughly the same value $`\gamma 0.7`$ for all continental records. Equation (1) can be used as a test bed for global climate models .

More recently, DFA was applied to study temperature correlations in the sea surface temperatures . It was found that the temperature autocorrelation function $`C(s)`$ again decayed by a power law, but with an exponent $`\gamma `$ close to 0.4, pointing towards a stronger persistence in the oceans than in the continents.

In this paper, we considerably extend our previous analysis to study systematically temperature records of 95 stations. Most of them are on the continents, and several of them are on islands. Our results are actually in line with both earlier papers and in agreement with conclusions drawn from independent type of analysis by several groups . We find that the continental records, including those on coastlines, show power-law persistence with $`\gamma `$ close to 0.7, while the island records show power-law correlations with $`\gamma `$ around 0.4. By comparing different orders of DFA that differ in the way trends are eliminated, we could also study the presence of trends in the records that lead to a warming of the atmosphere. We find that pronounced trends occur mainly at big cities and can be probably attributed to urban growth. Trends that cannot be attributed to urban growth occur in half of the island stations considered and on summit stations in the Alps. A majority of the stations showed no indications of trends.

The article is organized as follows. In Sec. II, we describe the detrending analysis used in this paper, the DFA. In Sec. III, we present the result of this analysis. Sec. IV concludes the paper with a discussion.

## II The methods of analysis

Consider a record $`T_i`$, where the index $`i`$ counts the days in the record, $`i=1,2,\mathrm{},N`$. The $`T_i`$ represent the maximum daily temperature, measured at a certain meteorological station. For eliminating the periodic seasonal trends, we concentrate on the departures of $`T_i`$, $`\mathrm{\Delta }T_i=T_i\overline{T}_i`$, from their mean daily value $`\overline{T}_i`$ for each calendar date $`i`$, say, 2nd of March, which has been obtained by averaging over all years in the record.

Quantitatively, correlations between two $`\mathrm{\Delta }T_i`$ values separated by $`n`$ days are defined by the (auto) correlation function

$$C(n)\mathrm{\Delta }T_i\mathrm{\Delta }T_{i+n}=\frac{1}{Nn}\underset{i=1}{\overset{Nn}{}}\mathrm{\Delta }T_i\mathrm{\Delta }T_{i+n}.$$
$`(2)`$
If $`\mathrm{\Delta }T_i`$ are uncorrelated, $`C(n)`$ is zero for $`n`$ positive. If correlations exist up to a certain number of days $`n_\times `$, the correlation function will be positive up to $`n_\times `$ and vanish above $`n_\times `$. A direct calculation of $`C(n)`$ is hindered by the level of noise present in the finite records, and by possible nonstationarities in the data.

To reduce the noise we do not calculate $`C(n)`$ directly, but instead study the “profile”

$$Y_m=\underset{i=1}{\overset{m}{}}\mathrm{\Delta }T_i.$$

We can consider the profile $`Y_m`$ as the position of a random walker on a linear chain after $`m`$ steps. The random walker starts at the origin and performs, in the $`i`$th step, a jump of length $`\mathrm{\Delta }T_i`$ to the right, if $`\mathrm{\Delta }T_i`$ is positive, and to the left, if $`\mathrm{\Delta }T_i`$ is negative. The fluctuations $`F^2(s)`$ of the profile, in a given time window of size $`s`$, are related to the correlation function $`C(s)`$. For the relevant case (1) of long-range power-law correlations, $`C(s)s^\gamma ,0<\gamma <1,`$ the mean-square fluctuations $`\overline{F^2(s)}`$, obtained by averaging over many time windows of size $`s`$ (see below) asymptotically increase by a power law :

$$\overline{F^2(s)}s^{2\alpha },\alpha =1\gamma /2.$$
$`(3)`$
For uncorrelated data (as well as for correlations decaying faster than $`1/s`$), we have $`\alpha =1/2`$.

For the analysis of the fluctuations, we employ a hierarchy of methods that differ in the way the fluctuations are measured and possible trends are eliminated (for a detailed description of the methods we refer to Ref. ).

(i) In the simplest type of fluctuation analysis (DFA0) (where trends are not going to be eliminated), we determine in each window the mean value of the profile. The variance of the profile from this constant value represents the square of the fluctuations in each window.

(ii) In the first order detrended fluctuation analysis (DFA1), we determine in each window the best linear fit of the profile. The variance of the profile from this straight line represents the square of the fluctuations in each window.

(iii) In general, in the $`n`$th order DFA (DFAn) we determine in each window the best $`n`$th order polynomial fit of the profile. The variance of the profile from these best $`n`$th-order polynomials represents the square of the fluctuations in each window.

By definition, DFA0 does not eliminate trends, while DFAn eliminates trends of order $`n`$ in the profile and $`n1`$ in the original time series. Thus, from the comparison of fluctuation functions $`F(s)`$ obtained from different methods one can learn about both, long-term correlations and the influence of trends.

The DFA method is analogous to wavelet techniques that also eliminate polynomial trends systematically. For a detailed review of the method, see Refs. . The conventional techniques such as the direct evaluation of $`C(n)`$, the rescaled range analysis (R/S) introduced by Hurst (for a review, see, e.g., Ref. ) or the power spectrum method  can only be applied on stationary records. In the presence of trends they may overestimate the long-term persistence exponent. The R/S method is somewhat similar to the DFA0 analysis.

## III Analysis of temperature records

Figure 1 shows the results of the DFA analysis of the daily temperatures (maximum or mean values) $`T_i`$ of the following weather stations (the length of the records is written within the parentheses): (a) Vienna (A, 125 yr), (b) Perm (RUS, 113 yr), (c) Charleston (USA, 127 yr), and (d) Pusan (KOR, 91 yr). Vienna and Perm have continental climate, while Charleston and Pusan are close to coastlines.

In the log-log plots the DFA1–3 curves are (except at small $`s`$ values) approximately straight lines. For both the stations inside the continents and along coastlines the slope is $`\alpha 0.65`$. There exists a natural crossover (above the DFA crossovers at very small times) that can be best estimated from DFA0 . As can be verified easily, the crossover occurs roughly at $`s_\times =10`$ days, which is the order of magnitude for a typical Grosswetterlage. Above $`s_\times `$, there exists long-range persistence expressed by the power-law decay of the correlation function with an exponent $`\gamma =22\alpha 0.7`$.

Figure 2 shows the results of the DFA analysis of the daily temperatures for two island stations: Wrangelja and Campbell Islands. Wrangelja Island is a large island between the East Siberian Sea and the Chukchi Sea. During the winter season, large parts of the water surrounding the island are usually frozen. Campbell Island is a small island belonging to New Zealand. Again, in the double logarithmic presentation, all DFA1–3 fluctuation functions are straight lines, but the slopes differ. While for Wrangelja the slope is 0.65, similar to the land stations shown before, the slope for Campbell Island is significantly larger, close to 0.8 (corresponding to $`\gamma =0.4`$).

It can be seen from Figs. 1 and 2 that sometimes the DFA0 curves have a larger slope than the DFA1–3 curves, and that usually the curves of DFA2 and DFA3 have the same slope for large $`s`$ values. The fact that the DFA0 curve has a higher exponent indicates the existence of trends by which the long-term correlations are masked. Calculations using DFA0 alone will yield a higher correlation exponent and thus lead to a spurious overestimation of the long-term persistence. The fact that the DFA2 and DFA3 curves show the same asymptotic behavior indicates that possible nonlinearities in the trends are not significant. Otherwise the DFA2 curve (where only linear trends are eliminated) would show an asymptotic behavior different from DFA3.

By comparing the DFA0 curves with the DFA2 curves, we can learn more about possible trends. Usually the effect of trends is seen as a crossover in the DFA0 curve. Below the crossover, the slopes of DFA0 and DFA2 are roughly the same, while above the crossover the DFA0 curve bends up. Large trends are characterized by a short crossovertime $`s_\times `$ and a large difference in the slopes between DFA0 and DFA2 (for a general discussion see Refs.  and ). A nice example for this represents Vienna, where the DFA0 curve shows a pronounced crossover at about 3 yr. Above this crossover, the DFA0 curve bends up considerably, with an effective slope close to 0.8. For Pusan, the trend is less pronounced, and for Perm, Charleston, and the two islands we do not see indications of trends.

To reveal the presence of long-term correlations and to point out possible trends, we have plotted in Fig. 3(a) the DFA0 curves and in Fig. 3(b) the DFA2 curves for 20 representative stations around the globe. For convenience, the fluctuation functions have been divided by $`s^{0.65}`$. We do not show results for those stations that were analyzed in Ref. . Figure 3(b) shows again that continental and coastline stations have roughly the same fluctuation exponent $`\alpha 0.65`$, while islands may also have higher exponents. It seems that stations at peaks of high mountains \[here we show Zugspitze (D, 98 yr, no. 19) and Wuxqiaoling (C, 40 yr, no. 20)\] have a slightly lower exponent.

From the 26 stations shown in Figs. 1–3, 8 show a larger exponent in the DFA0 treatment than in the DFA2 treatment. These stations are Thursday Island (no. 1 in Fig. 3), Koror Island (no. 2 in Fig. 3), as well as Vienna \[Fig. 1(a)\], Pusan \[Fig. 1(d)\], Hong Kong (no. 4 in Fig. 3), Jena (no. 10 in Fig. 3), Cheyenne (no. 17 in Fig. 3), and Zugspitze (no. 19 in Fig. 3). The other 18 stations do not show a difference in the exponents for DFA0 and DFA2, which suggests that the trends are either zero or too small to be detected by this sensitive method. We observe the largest trends for Hong Kong, Vienna, and Jena, where in all cases the crossover in the DFA0 curve is around 3 yr and the final slope is between 0.75 and 0.8. It is obvious that the greatest part of this warming is due to the urban growth of theses cities. Regarding the two islands, Koror shows a pronounced trend with a crossover time below 1 yr, while the trend we observe for Thursday Island is comparatively weak. It is not likely that the trends on the islands can be attributed to urban warming.

Figure 4 summarizes our results for all the stations analyzed. Fig. 4(a) shows the histogram for the values of the exponent $`\alpha `$ obtained by DFA0, while Fig. 4(b) shows the corresponding histogram obtained by DFA2. Both histograms are quite similar. For DFA2 the average exponent $`\alpha `$ is $`0.66\pm 0.06`$ and for DFA0 it is $`0.68\pm 0.07`$. The maxima become sharper when the islands are eliminated from the figures. The slight shift towards larger $`\alpha `$ values in DFA0 is due to trends. The magnitude of the trends can be roughly characterized by the difference $`\delta \alpha `$ of the slopes of DFA0 and DFA2. We found that 7 of the 15 island stations and 54 of the 80 continental stations showed no significant trend, with $`\delta \alpha 0.02`$. We observed a small trend, with $`0.03\delta \alpha 0.05`$, for 3 island and 9 continental stations. A pronounced trend, with $`\delta \alpha 0.06`$, was found for 5 island and 13 continental stations. Among these 13 stations are Hong Kong, Bordeaux, Prague, Seoul, Sydney, Urumchi, Swerdlowsk, and Vienna, where a large part of the warming can be attributed to the urban growth of the cities in the last century. Two of these stations \[Säntis (CH) and Sonnblick (A)\] are on top of high mountains.

Since the island stations have exponent $`\alpha `$ larger than the continental stations, it is likely that the long-term persistence originates by the coupling of the atmosphere to the oceans. Thus one might expect that for island stations $`\alpha `$ will increase with the distance to the continents, and for continental stations $`\alpha `$ will decrease with the distance to the coastline. To test if the exponent $`\alpha `$ depends on the distance $`d`$ to the continental coastlines, we plotted in Fig. 5 $`\alpha `$ as a function of $`d`$ for both islands and continental stations. It is remarkable that islands far away from the continents do not show a larger exponent than islands close to the coastlines, and inner-continental stations far from the ocean do not show smaller exponents than coastline stations. This second result is in disagreement with a recent claim that $`\alpha =0.5`$ for inner-continental stations far away from the oceans .

## IV Discussion

In this paper, we have used a hierarchy of detrending analysis methods (DFA0–DFA3) to study long temperature records around the globe. We concentrated mainly on those areas on the globe (North America, Europe, Asia and Australia) where long records are available. The main results of the study are the following.

(i) The temperature persistence decays, after a crossover time that is typically of the order of the duration of a Grosswetterlage, by a power law, with an exponent $`\alpha `$ that has a very narrow distribution for continental stations. The mean value of the exponent is close to 0.65, in agreement with earlier calculations based on different methods .

(ii) On islands, the exponent shows a broader distribution, varying from 0.65 to 0.85, with an average value close to 0.8. This finding is in qualitative agreement with the results of a recent analysis of sea surface temperature records, where also long-term persistence with an average exponent close to 0.8 has been found . Since the oceans cover more than 2/3 of the globe, one may expect that also the mean global temperature is characterized by long-term persistence, with an exponent close to 0.8.

(iii) In the vast majority of stations we did not see indications for a global warming of the atmosphere. Exceptions are mountain stations in the Alps \[Zugspitze (D), Säntis (CH), and Sonnblick (A)\], where urban warming can be excluded. Also, in half of the islands we studied, we found pronounced trends that most probably cannot be attributed to urban warming. Most of the continental stations where we observed significant trends are large cities where probably the fast urban growth in the last century gave rise to temperature increases.

When analyzing warming phenomena in the atmosphere, it is essential to employ methods that can distinguish, in a systematic way, between trends and long-term correlations in contradistinction to a number of conventional schemes that have been applied in the past. These schemes run the risk of mixing up the correlatedness of natural climate system variability with entire regime shifts enforced by anthropogenic interference through greenhouse gas emissions. The fact that we found it difficult to discern warming trends at many stations that are not located in rapidly developing urban areas may indicate that the actual increase in global temperature caused by anthropogenic perturbation is less pronounced than estimated in the last IPCC (Intergovernmental Panel for Climate Change) report .

## Acknowledgments

We are grateful to Professor S. Brenner for very useful discussions. We would like to acknowledge financial support by the Deutsche Forschungsgemeinschaft and the Israel Science Foundation.


format: text

text:
# Thermal stability of coupled ferromagnetic and superparamagnetic particles

## I Introduction

Thermal stability is of great importance in ultra-high-density magnetic recording. The principal element of magnetic memory is a single-domain ferromagnetic particle. The magnetic moment of this particle is agitated by thermal fluctuations and can reverse to the opposite direction. In the absence of inter-particle interactions the magnetization reversal time is determined by the Néel-Arrhenius exponent with a definite energy barrier (e.g., ,,). The problem for the mean magnetization reversal time has been rigorously solved in the framework of a new theoretical approach, in which the coherent rotation of the magnetization is analyzed in terms of the motion of a damped nonlinear oscillator subject to random Gaussian forces (thermal agitation) ,, ,,.

Recently it has been experimentally demonstrated that the thermal stability of magnetic memory can increase if the ferromagnetic particles from the recording layer are antiferromagnetically coupled to an auxiliary layer of smaller ferromagnetic particles (AFC media) ,. Simple theoretical explanations of this phenomenon are given, for example, in Refs.,. The basic model consists of two single-domain hard particles (recording and stabilizing) coupled by an antiferromagnetic exchange interaction.

Here we analyze the possibilities of the two-particle model to enhance thermal stability of recording media. We consider the thermally agitated magnetization reversal, when a hard ferromagnetic particle interacts with a soft ferromagnetic (superparamagnetic) particle. We use the new theoretical approach for coupled damped nonlinear oscillators subject to random Gaussian forces introduced by thermal agitation.

## II Model

Let us consider a ferromagnetic particle with magnetization $`𝐌_1`$($`|𝐌_1|=M_{s,1}`$) and volume $`V_1`$ and a superparamagnetic particle with $`𝐌_2`$ ($`|𝐌_2|=M_{s,2}`$) and $`V_2`$, respectively (Fig.1). To describe the coherent motion of the particle magnetizations we introduce the classical spins $`𝐒_j=𝐌_jV_j/\mathrm{}\gamma `$ ($`j=1,2`$), where $`\gamma `$ is the gyromagnetic ratio and $`\mathrm{}`$ is Planck’s constant.

The energy of the combined system includes uniaxial anisotropy of the ferromagnetic particle along $`z`$ axis, Zeeman energy, exchange and magnetostatic interparticle interactions:

$``$ $`=`$ $`K_uV_1\left({\displaystyle \frac{S_1^z}{S_1}}\right)^2\mathrm{}\gamma 𝐇_0(𝐒_1+𝐒_2)+J𝐒_1𝐒_2`$ (2)
$`+(\mathrm{}\gamma )^2\left[{\displaystyle \frac{𝐒_1𝐒_2}{r^3}}3{\displaystyle \frac{(𝐒_1𝐫)(𝐒_2𝐫)}{r^5}}\right].`$
Here $`K_u`$ is the uniaxial anisotropy constant, $`𝐇_0`$ is the external magnetic field and $`J`$ is the exchange integral ($`J>0`$ corresponds to antiferromagnetic and $`J<0`$ to the ferromagnetic coupling, respectively). The last term in (2) describes the magnetostatic interactions (dipole-dipole approximation), where $`𝐫`$ is the vector between particle centers.

For simplicity we shall confine our analysis to two particle configurations applicable for the longitudinal (Fig.1a) and perpendicular (Fig.1b) recording media with magnetic fields parallel to each uniaxial anisotropy axis. Thus Eq.(2) becomes

$``$ $`=`$ $`K_uV_1\left({\displaystyle \frac{S_1^z}{S_1}}\right)^2\mathrm{}\gamma H_0(S_1^z+S_2^z)`$ (4)
$`+\mathrm{}\gamma (H_{int}S_1^zS_2^z+H_{,x}S_1^xS_2^x+H_{,y}S_1^yS_2^y).`$
Here $`H_{,y}=J/\mathrm{}\gamma +\mathrm{}\gamma /r^3`$. For configurations shown in Fig.1 the interaction fields $`H_{int}`$ and $`H_{,x}`$ are defined by

$`H_{int}^{(a)}`$ $`=`$ $`{\displaystyle \frac{J}{\mathrm{}\gamma }}+{\displaystyle \frac{\mathrm{}\gamma }{r^3}},`$ (5)
$`H_{,x}^{(a)}`$ $`=`$ $`{\displaystyle \frac{J}{\mathrm{}\gamma }}2{\displaystyle \frac{\mathrm{}\gamma }{r^3}}`$ (6)
for Fig.1a and

$`H_{int}^{(b)}`$ $`=`$ $`{\displaystyle \frac{J}{\mathrm{}\gamma }}2{\displaystyle \frac{\mathrm{}\gamma }{r^3}},`$ (7)
$`H_{,x}^{(b)}`$ $`=`$ $`{\displaystyle \frac{J}{\mathrm{}\gamma }}+{\displaystyle \frac{\mathrm{}\gamma }{r^3}}`$ (8)
for Fig.1b.

## III Dynamic equations

In order to describe spin (magnetization) dynamics, we introduce complex variables $`a_j^{}`$ and $`a_j`$ (classical analog of creation and annihilation operators) of nonlinear oscillators defined by the Holstein-Primakoff transformation :

$`S_j^z`$ $`=`$ $`S_jN_j,N_j=a_j^{}a_j,`$ (9)
$`S_j^x`$ $`=`$ $`{\displaystyle \frac{a_j+a_j^{}}{2}}\sqrt{2S_jN_j},`$ (10)
$`S_j^y`$ $`=`$ $`{\displaystyle \frac{a_ja_j^{}}{2i}}\sqrt{2S_jN_j}.`$ (11)
This transformation gives a convenient representation of spins as oscillators in the vicinity of equilibrium $`𝐒_1=(0,0,S_1)`$ and $`𝐒_2=(0,0,S_2)`$.

The dynamic equations for $`a_1`$ and $`a_2`$, defined by $`da_j/dt=i(/\mathrm{})/a_j^{}`$, can be written as

$`{\displaystyle \frac{d}{dt}}a_1`$ $`=`$ $`i\gamma H_{eff,1}a_1i\gamma F_{12},`$ (12)
$`{\displaystyle \frac{d}{dt}}a_2`$ $`=`$ $`i\gamma H_{eff,2}a_2i\gamma F_{21}.`$ (13)
Here

$`H_{eff,1}`$ $`=`$ $`H_0+H_KS_1^zH_{int}S_2^z,`$ (14)
$`H_{eff,2}`$ $`=`$ $`H_0H_{int}S_1^z`$ (15)
are the effective fields, $`H_K=2K_u/M_{s,1}`$ is the anisotropy field and

$$F_{mn}=H_{,x}\frac{S_n^x}{a_n^{}}S_m^x+H_{,y}\frac{S_n^y}{a_n^{}}S_m^y.$$
(16)
The principal oscillatory motion of $`a_1(t)`$ and $`a_2(t)`$ are determined by their effective fields: $`a_1\mathrm{exp}(i\gamma H_{eff,1}t)`$ and $`a_2\mathrm{exp}(i\gamma H_{eff,2}t)`$. Because, in general, $`H_{eff,1}H_{eff,2}`$, the fast oscillating terms $`F_{12}`$ and $`F_{21}`$ in (12) vanish. The magnetization fluctuations of the superparamagnetic particle occurs at a much faster rate $`f_2`$ than the reversal rate of the ferromagnetic particle. $`f_2`$ is an attempt frequency of the second particle (a negligibly small energy barrier is assumed). Thus, the longitudinal component $`S_1^z`$ is the slowest variable in the system, and the analysis of the coupled system can be reduced to the stochastic dynamics of the first particle in an averaged field from the second particle. This field appears as a feedback of the superparamagnetic particle on the instantaneous ordered state $`S_1^z`$ of the first particle. The corresponding effective energy is equal to

$$_2\mathrm{}\gamma H_{eff,2}S_2^z.$$
(17)
In this mean field approximation (17), it follows that for times $``$ $`f_2^1`$ an averaged thermal spin polarization of the second (superparamagnetic) particle (e.g. ) is:

$$S_2^z=S_2B_{S_2}\left(S_2\frac{\mathrm{}\gamma H_{eff,2}}{k_BT}\right),$$
(18)
where

$$B_S(Sx)=1\frac{1}{S}\left(\frac{1}{\mathrm{exp}(x)1}\frac{2S+1}{\mathrm{exp}[(2S+1)x]1}\right)$$
(19)
is the Brillouin function. This averaged spin polarization $`S_2^z`$ induces an effective field (or, response field) back to the first particle. As a result the effective energy of the first particle can be written as

$$_1\frac{\mathrm{}\gamma H_K}{2S_1}\left(S_1^z\right)^2\mathrm{}\gamma \left(H_0H_{int}S_2^z\right)S_1^z.$$
(20)
From (5) and (7) one can see that $`H_{int}`$ is enhanced for antiferromagnetic ($`J>0`$) exchange in the longitudinal case (Fig.1a) and for ferromagnetic exchange in the perpendicular case (Fig.1b).

## IV Approximate solution

Let us consider a small polarization of the superparamagnetic particle. In this case Eq.(18) is simplified to

$$S_2^z\frac{S_2(S_2+1)}{3}\frac{\mathrm{}\gamma H_{eff,2}}{k_BT}$$
(21)
and the effective energy (20) becomes

$$_1\frac{\mathrm{}\gamma \stackrel{~}{H}_K}{2S_1}\left(S_1^z\right)^2\mathrm{}\gamma H_{eff,1}S_1^z.$$
(22)
Here

$$\stackrel{~}{H}_K=H_K\left(1+2S_1\zeta \frac{H_{int}}{H_K}\right)$$
(23)
is the effective anisotropy field and

$$H_{eff,1}=H_0\left(1\zeta \right)$$
(24)
is the effective external field on the first particle,

$$\zeta =\frac{S_2(S_2+1)}{3}\frac{\mathrm{}\gamma H_{int}}{k_BT}.$$
(25)
From (23) and (25) we see that the effective anisotropy field $`\stackrel{~}{H}_K`$ always increases ($`H_{int}\zeta H_{int}^2`$) for any $`H_{int}0`$. On the contrary, the effective external magnetic field $`H_{eff,1}`$ (24) depends on the sign of $`H_{int}`$: it can be both enhanced ($`\zeta H_{int}<0`$) and suppressed ($`\zeta >0`$). The feedback effect from the superparamagnetic particle quadratically depends on the saturation magnetization and volume of this particle: $`\zeta S_2^2M_2^2V_2^2`$. This quadratic dependence also means that the corresponding feedback effect of two superparamagnetic particles of volume $`V_2/2`$ is at least two times smaller. Thus it is more effective to use a single larger superparamagnetic particle than a collection of smaller ones.

The energy (22) is an exact analog to the case of random motion of one single-domain particle when the external magnetic field is parallel to the anisotropy axis. Thus, we can use exact analytic formulas ,, for the mean first passage time obtained for this case. It should be noted that this approach makes it possible to describe experimental data for the reversal field versus pulse time over 13 orders from nanoseconds to hours , including both purely dynamic and thermally assisted reversal.

## V Results and discussion

Let us estimate Eqs. (23) and (24). For the ferromagnetic particle we take the following parameters at room temperature: $`K_uV_1/k_BT=60`$, $`H_K=5`$ kOe and $`M_{s,1}=150`$ emu/cc. This means $`V_16.6210^{18}`$ cm3 and $`S_15.3610^4`$. Taking $`S_2700`$ and $`H_{int}=0.3`$ Oe, one obtains: $`\zeta 0.022`$ and $`2S_1\zeta H_{int}/H_K0.14`$. Thus, we see that the effect of the superparamagnetic particle is relatively small in the effective magnetic field $`H_{eff,1}0.978H_0`$ and larger in the effective anisotropy field $`\stackrel{~}{H}_K=1.14H_K`$. For the opposite sign of interaction (ferromagnetic exchange between particles) $`H_{int}=0.3`$ Oe and with the same parameters, one has $`H_{eff,1}1.022H_0`$ and $`\stackrel{~}{H}_K=1.14H_K`$.

In Fig.2 reversal field versus pulse time is shown for the above two cases. We also show the corresponding curve for one ferromagnetic particle ($`H_{int}=0`$) and Sharrock’s (Néel-Arrhenius) dependence. The thick solid line and broad dashed curve represent antiferromagnetic and ferromagnetic coupling, respectively. The slight difference is due to the larger effective field and hence reduced stability of the ferromagnetic case. In the thermal region ($`H<H_K`$) there is more curvature than in the Néel-Arrhenius approximation. This is because the nonlinear oscillator model contains a field dependent attempt frequency. Note that formulas (22)-(25) are valid for long time scale $``$ $`f_2^110^{10}10^9`$ sec. In the dynamic reversal region particles seems to reverse independently.

## Acknowledgment

This work was partly supported by matching funds from the Center for Magnetic Recording Research at the University of California - San Diego and CMRR incorporated sponsor accounts.

Figure captions

Fig.1. Particle configurations for a) longitudinal and b) perpendicular recording. The larger volume represents the ferromagnetic particle.

Fig.2. Reversal field versus pulse time.


text:
# Magnetic phase diagram of Ca1-xMnxO

## Abstract

Alternating current susceptibility and direct current magnetization have been studied for polycrystalline Ca1-xMnxO. On increasing the Mn content, magnetic ordering changes from spin-glass behavior for $`0.25x0.4`$ to antiferromagnetic order. The paramagnetic/antiferromagnetic transition is of second order for $`0.5x0.65`$ and of first order for $`x0.7`$. For low Mn concentrations, the high-temperature alternating current susceptibility can be described by a diluted Heisenberg magnet model developed for diluted magnetic semiconductors.

Theoretical predictions of room temperature ferromagnetism in Mn-containing diluted magnetic semiconductors Dietl00  recently brought wide attention to this class of materials. According to these calculations, $`p`$-type Zn1-xMnxO is a promising candidate for a room temperature ferromagnet. Pulsed-laser deposited Zn1-xMnxO thin films without intentional carrier doping show spin-glass behavior.Fukumura01  According to this latter study, Mn can be dissolved in the ZnO matrix to over 35%. We have shown in a previous paper that the solubility of Mn in the zinc-blende structure of ZnO is less than 15% for polycrystalline samples under various oxygen pressure conditions. Samples obtained in air or argon are paramagnetic, while the high-pressure oxygen annealing induces spin-glass-like behavior by precipitation of ZnMnO3 in the paramagnetic matrix.Kolesnik02

In this study we investigate polycrystalline Ca1-xMnxO $`(x=0.1251)`$. Both end-member compounds, CaO and MnO, adopt the rock-salt structure and the complete range of solid solution can be achieved. This makes it possible to study a continuous evolution of magnetic properties from diamagnetic insulating/semiconducting CaO (energy gap $`E_g=7.09`$ eV)Whited69  to antiferromagnetic insulating MnO ($`T_N120`$ K, $`E_g=4`$ eV), whose antiferromagnetic properties were already observed seventy years ago.Tyler33  Due to similarity of the ionic radii (0.8 Å for Mn2+ and 0.99 Å for Ca2+) Mn2+ ions should occupy a substitutional site in CaO in contrast to BaO (ionic radius of Ba2+ is equal to 1.34 Å), where Mn2+ ions locate off-center.Winsum78  The CaO-MnO solid solution has been shown to form in the H2/He atmosphere above 1173 K from carbonate precursors.Poeppelmeier86  Thermodynamic functions of mixing have been determined at high temperatures for Ca1-xMnxO.Rog92  To our knowledge, the magnetic properties of this solid solution have not yet been reported. MnO undergoes a first-order transition to the antiferromagnetic state at $`T_N=115120`$ K Seino73 ; Srinivasan83  This transition is accompanied by a trigonal deformationMorosin70  to a slightly rhombohedral structure, compressed along the  direction.Bloch73  The antiferromagnetic order, determined by neutron diffraction, is of type-II (AFII) in the fcc lattice and consists of (111) sheets of ferromagnetically ordered Mn2+ ions, which are antiferromagnetically coupled to the neighboring sheets along the  direction.Roth58  Magnetic properties of MnO have been explained in terms of two superexchange constants $`J_1=10`$ K and $`J_2=11`$ K for the nearest-neighbor (NN) and next-nearest-neighbor (NNN) interactions, respectively.Lines65  Recent local-density approximation (LDA) calculationsPask01  suggest relatively stronger NNN interactions $`J_1=9.8`$ K and $`J_2=24.5`$ K. The small uniaxial stress along the  axis of the MnO crystal results in the clearly first-order transition and the increasing $`T_N`$.Bloch73  This result has been ascribed to the strain-related suppression of T domains (which correspond to antiferromagnetic stacking along the four  directions within the crystal). Due to the presence of these T domains the first-order character of the antiferromagnetic/paramagnetic transition may be obscured. Recent heat capacity studies showed a continuous magnetic transition in small, stress-free, MnO crystals.Woodfield99  The application of large  stress ($`1.2<\tau 5.5`$ kbar) leads to a continuous transition, again. Bloch75  Solid solution Mn1-xNixO shows a linear increase of $`T_N`$ upon increasing Ni concentration with preservation of the AFII structure.Cheetham83  Both end members of this series, MnO and NiO ($`T_N525`$ K) show the same type of rhombohedral crystallographic distortions.Roth58

The mixture of CaCO3 and MnO2 was calcined in air several times at $`T=900,1000,1100^{}`$C with intermediate grindings. The final synthesis stage was performed in a hydrogen atmosphere at $`T=1050^{}`$C. The alternating current (ac) susceptibility and direct current (dc) magnetization were measured using a Physical Property Measurement System (Quantum Design). The ac susceptibility was measured in the range 2.5-395 K, in the excitation field $`H_{ac}=14`$ Oe at the frequency $`\omega =1`$ kHz. The “zero-field-cooled” ($`M_{ZFC}`$) and “field-cooled” ($`M_{FC}`$) magnetizations were measured in a magnetic field of 1 kOe. $`M_{ZFC}`$ was measured on warming after cooling in a zero magnetic field and switching the magnetic field on at $`T=5`$ K. $`M_{FC}`$ was subsequently measured on cooling in the magnetic field.

X-ray diffraction spectra have been collected using a Rigaku X-ray diffractometer. The sharp X-ray diffraction patterns, presented in Fig. 1 for several samples, indicate good quality, homogeneous material.

The single-phase rock-salt-type structure is revealed for all the Ca1-xMnxO samples. Inset to Fig. 1 shows the lattice parameter, $`a`$, determined from the X-ray diffraction patterns using the GSAS refinement software. The lattice parameter, $`a`$, is linearly dependent on the Mn content $`x`$ with a very good accuracy. The linear fit to the experimental data gives a dependence $`a(x)=4.812(1)0.369(2)x`$ (Å). For CaO, the extrapolated $`a`$ is equal to 4.812(1) Å, which is in agreement with the literature value $`a=4.81`$ Å.McMurdie86  For our MnO sample, $`a=4.443(3)`$ Å is very close to the literature value $`a=4.4424.446`$ Å.Barrett64  A very good agreement between the measured lattice parameters and the linear dependence for all Mn concentrations suggests that the actual compositions and the nominal ones are the same.

Three different kinds of magnetic transitions observed for Ca1-xMnxO are illustrated in

Fig. 2. For $`x=0.250.4`$, the spin-glass behavior is observed at low temperatures. A difference between the “zero-field-cooled” ($`M_{ZFC}`$) and “field-cooled” ($`M_{FC}`$) magnetizations, which is characteristic of spin-glass behavior, is shown in Fig. 2(a) for $`x=0.4`$. The presence of the thermoremanent ($`M_{trm}`$) magnetization at low temperatures supports this observation. The inset to Fig. 2(a) shows a cusp in the ac susceptibility. The temperature, at which the cusp is maximum defines the spin-glass freezing temperature, $`T_f`$. Similar cusp was observed for the $`x=0.333`$ sample and a kink on the ac susceptibility for $`x=0.25`$. The same type of spin-glass behavior has been previously observed for Mn-containing diluted magnetic semiconductors (e.g. Zn1-xMnxSe Twardowski87  and Cd1-xMnxTe Denissen87 ), where a cusp was observed for high Mn concentrations and a kink for low Mn concentrations.

For $`x=0.50.65`$, $`M_{ZFC}`$ and $`M_{FC}`$ curves trace each other and no thermoremanent magnetization is observed. This behavior, presented in Fig. 2(b) indicates a continuous (second-order) paramagnetic/antiferromagnetic transition. For higher Mn concentrations, $`x0.7`$ including MnO, the paramagnetic/antiferromagnetic transition is accompanied by a hysteresis, clearly observed both in the ac susceptibility and dc magnetization (see: Fig. 2(c)). This hysteresis is a sign of the discontinuous first-order transition.

The observed magnetic transition temperatures are collected in Fig. 3.

The $`x=0.125`$ sample is paramagnetic down to $`T=2.5`$ K. The $`x=0.250.4`$ samples show spin-glass behavior at increasing temperatures with increasing Mn content. The antiferromagnetic transition temperatures for $`x0.5`$ are defined as the temperatures for which the temperature dependence of the ac susceptibility has the maximum slope, i.e. slightly below the observed maximum of $`\chi (T)`$. For the first-order transitions, ($`x0.7`$), we used the transition temperatures, which were determined on warming ($`M_{ZFC}`$). The transition temperatures in all the Mn concentration regions seem to continuously change with $`x`$, but no common scaling has been determined. The spin-glass behavior at low Mn concentrations is quantitatively similar to the behavior of wurtzite and zinc-blende diluted magnetic semiconductors. The temperature dependence of the ac susceptibility for $`x0.5`$ is almost identical with that of the dc magnetization. The width of the thermal hysteresis $`\mathrm{\Delta }T`$ in the ac susceptibility at the antiferromagnetic transition, which is shown in the inset to Fig. 3 starts increasing with Mn concentration for $`x2/3`$. This indicates a crossover between the second-order and first-order nature of the transition.

The high-temperature ac susceptibility of Ca1-xMnxO generally shows Curie-Weiss behavior. Within the framework of the diluted Heisenberg antiferromagnet theory,Spalek86  the high-temperature ac mass susceptibility can be described by the formula

$$\chi =\frac{C_M(x)}{\mu (x)(T\mathrm{\Theta }(x))},$$
(1)
where $`\chi `$ is the ac susceptibility after subtraction of the diamagnetic contribution of the host material CaO (equal to $`0.2710^6`$ emu/g) and $`C_M(x)`$ is the molar Curie constant defined as

$$C_M(x)=x\frac{N(g\mu _B)^2S(S+1)\mu (x)}{3k_B\rho (x)},$$
(2)
$`N`$ is the number of cations per unit volume, $`S`$ is the effective spin of Mn2+ ion, $`g=2`$ is the gyromagnetic factor, $`\mu (x)`$ is the molar mass, $`\rho (x)`$ is the mass density, $`\mathrm{\Theta }(x)=\mathrm{\Theta }_0x`$ is the Curie-Weiss temperature

$$\mathrm{\Theta }_0(x)=\frac{2}{3}xS(S+1)\underset{P}{}J_Pz_P/k_B\mathrm{\Theta }_0x,$$
(3)
where $`J_P`$ is the exchange integral for a pair of $`P`$th Mn neighbors and $`z_P`$ is the number of neighbors of $`P`$th order. In the diluted magnetic regime the dominating interactions are between the nearest neighbors. Hence, we may limit the expansion in Eq. (3) to the first order. The $`\mathrm{\Theta }_0`$ constant is then related to the exchange integral between the nearest Mn neighbors $`J_1`$,

$$\frac{2J_1}{k_B}=\frac{3\mathrm{\Theta }_0}{zS(S+1)},$$
(4)
where $`z=12`$ is the number of nearest neighbors in the rock-salt structure of Ca1-xMnxO.

The inverse ac susceptibility for several Ca1-xMnxO samples is shown in Fig. 4.

For low Mn concentrations, $`x0.333`$, the inverse ac susceptibility is linear with temperature in a wide range $`T=50390`$ K. Deviations from linearity (characteristic to diluted antiferromagnets and related to the interactions between next nearest neighboring Mn2+ ionsSpalek86 ) are observed at lower temperatures. For higher Mn concentrations, the slope of the inverse susceptibility changes with temperature. We have analyzed the linear fits to the $`\chi ^1(T)`$ data in three temperature ranges: 100-200 K, 200-300 K, and 300-390 K.

The data points in Fig. 5 represent the parameters determined from the fit of Eq. (1) to the inverse susceptibility data in the temperature range 200-300 K.

The values determined from the fits in the other two temperature regions are plotted as error bars. For a given $`x`$, the fit in the lower temperature range (100-200 K) gives higher absolute $`\mathrm{\Theta }(x)`$ and larger $`C_M(x)`$ than the fit in the higher temperature range (300-390 K).

For diluted Heisenberg antiferromagnets it is expected that the material parameters $`C_M(x)`$ and $`\mathrm{\Theta }(x)`$ are proportional to $`x`$. We have made linear fits to the first three data points (i.e. $`0.125x0.333`$, assuming that $`C_M(x=0)=0`$ and $`\mathrm{\Theta }(x=0)=0`$). From these fits, we have obtained the extrapolated values $`C_M(x=1)=4.53(4)`$ emu K/mol and $`\mathrm{\Theta }_0=303(11)`$ K. The value of $`C_M`$ is close to those determined for Mn-containing II-VI diluted magnetic semiconductors.Spalek86  From Eq. (2), calculated $`S=2.55(5)`$, which is in perfect agreement with the effective spin of Mn2+ ion $`S=5/2`$. The extrapolated $`\mathrm{\Theta }_0`$ is lower (in the sense of the absolute value) than the values found for diluted magnetic semiconductors ($`\mathrm{\Theta }_0=470`$ to $`831`$ K).Spalek86  We observe significant deviations from the linear behavior of $`\mathrm{\Theta }(x)`$ for $`x0.5`$, where the antiferromagnetic order can be observed and the diluted magnet model is no longer valid. It is worth noting that for higher Mn concentrations the Curie-Weiss temperatures become comparable to or higher than our experimental temperature range, which increases the error of the determined values. However, we suggest that these deviations of $`\mathrm{\Theta }(x)`$ from the linear behavior can be explained by the increasing contribution of NNN interactions between highly concentrated Mn2+ ions.

By including the NNN interactions we obtain the following modification of Eq. (3)

$$z\frac{2J_1}{k_B}+z_2\frac{2J_2}{k_B}=\frac{3\mathrm{\Theta }(x)}{S(S+1)},$$
(5)
where $`z_2=6`$ is the number of next-nearest neighbors. For MnO, for example, $`\mathrm{\Theta }=536`$ K. By using $`2J_1/k_B=8.7(3)`$ K calculated from the extrapolated $`\mathrm{\Theta }_0=303`$ K according to Eq. (4), we obtain $`2J_2/k_B13.2`$ K for MnO, which gives us a qualitative agreement with the LDA calculationsPask01 , namely $`J_2>J_1`$. However, the actual values of our exchange constants are significantly smaller than those found in the literature (see: references to Ref. Pask01 ). Note that M. E. Lines and E. D. JonesLines65  used an analogical formula to our Eq. (5) but without the factor 2. From their $`\mathrm{\Theta }=540`$ K, this gave as a result $`(2J_2+J_1)/k_B=30.9`$ K, which is twice as higher as our value.

In summary, we have determined the magnetic phase diagram of Mn- doped polycrystalline CaO samples. On increasing the Mn content, magnetic ordering changes from the spin-glass behavior for $`0.25x0.4`$ to the antiferromagnetic order. The paramagnetic/antiferromagnetic transition is second order for $`0.5x0.65`$ and first order for $`x0.7`$. For low Mn contents, the ac susceptibility can be described by a diluted Heisenberg magnet model developed for diluted magnetic semiconductors.

This work was supported by the DARPA/ONR and the State of Illinois under HECA.


text:
# A nonlinear dynamical model of human gait

## I Introduction

Walking is a complex process which we have only recently begun to understand through the application of nonlinear data processing techniques hausdorff1 ; hausdorff99 ; west1 ; west2  to study interval data. It has been known for over a century that there is variation of 3-4% in the stride interval of humans during walking vierordt , and only in the last decade did Hausdorff et al. hausdorff1  demonstrate that the stride-interval time series exhibits long-time correlation, suggesting that the phenomenon of walking is a self-similar, fractal activity. Subsequent studies by West and Griffin west1 ; west2  supported the conclusion that the human gait time series is fractal. However, more recently it was determined that these time series, rather than being monofractal, are weakly multifractal griffin ; scafetta1 .

Human locomotion is known to be a voluntary process, but it is also regulated through a network of neurons called a Central Pattern Generator (CPG) collins1 , capable of producing a syncopated output. The early nonlinear dynamical models of CPGs for gait assumed that a single nonlinear oscillator be used for each limb participating in the locomotion process collins2 . Therefore a quadraped required the coupling of four nonlinear oscillators to determine the correct phase relations among the four legs in order to distinguish between various modes of locomotion, that is, walking, trotting, cantering and galloping. More recent dynamical models, using the property of synchronization of nonlinear dynamical systems, allow for neurons within an assembly to become enslaved to a single rhythmic muscular activity. Thus, rather than having a separate nonlinear oscillator for each limb, it is possible to have a single CPG to determine how we walk.

The model that we present here, the super CPG (SCPG), assumes that the central nervous system is coupled to the motocontrol system, and together they control the locomotion of the gait cycle. We stress that it is the period of the gait cycle that is ultimately measured in these stride interval experiments, and not the neural firing activity. The dynamics of human gate may also be voluntarily forced, for example, by following the frequency of a metronome. We model the complex gait system by assuming that the amplitude of the impulses of the correlated firing neural centers regulate only the unperturbed inner frequency of a nonlinear forced Van der Pol oscillator holmes  that mimics the gait cycle. The stride interval is assumed to coincide with the actual period of the Van der Pol oscillator. In this way the gait frequency may differ slightly from the potential frequency induced by the neural firing activity. In fact, the chaotic behavior of nonlinear oscillators, like the Van der Pol oscillator, allows a more complex behavior that may be controlled also by an constraint that forces the oscillator to follow a particular fixed frequency.

The SCPG model is tested on a data set available in the public domain archives Physionet physionet . These data were originally collected and used by Hausdorff et al. Hausdorff4  to determine the dependence of the fractal dimension of the time series on changes of the average rate of walking. These data contain the stride interval time series for 10 healthy young men walking at a slow, normal and fast pace, for a period of one hour. The same individuals are then requested to walk at a pace determined by a metronome set at the average slow, normal and fast paces for 30 minutes to generate a second data set.

By estimating the Hölder exponents and their spectra using wavelet transform struzik , it was shown that the stride interval time series is weakly multifractal with a main fractality close to that of the 1/f-noise. The time series is sometimes non-stationary and its fractal variability changes in the different gait mode regimes scafetta1 . In particular, the persistence, as well as the multifractality of the stride interval time series, tends to increase for both slow and fast pace, above that of the normal paces. Moreover, if the pace is constrained by a metronome, the stochastic properties of the stride interval time series change significantly, from persistent to antipersistent fluctuations, but, in general, in each case there is a reduction in the long-term memory and an increase in randomness.

In Sec. 2 we give a short introduction to the phenomenon of locomotion, the traditional methods for modeling using the CPG, and review the data processing used to establish the fractal behavior of stride interval time series. Sec. 3 reviews the stochastic properties of the normal and metronomic gait under different various pace velocities, slow, normal and fast. In Sec. 4 we present the mathematical details of the SCPG model. In Sec. 5 we compare the results of computation using the SCPG model with the phenomenological data. Finally, in Sec. 6 we draw some conclusions.

## II Central pattern generator and locomotion

Walking consists of a sequence of steps. These steps may be partitioned into two phases: a stance phase and a swing phase. The stance phase is initiated when a foot strikes the ground and ends when it is lifted. The swing phase is initiated when the foot is lifted and ends when it strikes the ground again. The time to complete each phase varies with the stepping speed. A stride interval is the length of time from the start of one stance phase to the start of the next stance phase.

Traditionally the legged locomotion of animals is understood through the use of a CPG, an intraspinal network of neurons capable of producing a syncopated output collins1 ; crago . The implicit assumption in such an interpretation is that a given limb moves in direct proportion to the voltage generated in a specific part of the CPG. Experiments establishing the existence of a CPG have been done on animals with spinal cord transections. Walking, for example, in a mesencephalic cat, a cat with its brain stem sectioned rostral to the superior colliculus, is very close to normal, on a flat, horizontal surface, when a section of the midbrain is electrically stimulated. Stepping continues as long as a train of electrical pulses is used to drive the stepping. This is not a simple linear response process because changing the frequency of the driver has little effect of the walking cycle mann . However, since the frequency of the stepping increases in proportion to the amplitude of the stimulation, we can conclude that the variation in the stride interval of humans is related to the fluctuation of the amplitude of the impulses of the firing neural centers.

As Collins and Richmond collins1  point out, in spite of the studies establishing the existence of a CPG in the central nervous system of quadrupeds, such direct evidence does not exist for a vertebrate CPG for biped locomotion. Consequently, these and other authors have turned to the construction of models, based on the coupling of nonlinear oscillators, the hard-wired CPG, to establish that the mathematical models are sufficiently robust to mimic the locomotion characteristics observed in the movements of segmented bipeds cohen , as well as in quadrupeds collins2 . These characteristics, such as the switching among multiple gait patterns, is shown to not depend on the detailed dynamics of the constituent nonlinear oscillators, nor on their inter-oscillator coupling strengths collins1 .

As we mentioned in the Introduction, it has been known for over a century that there is a variation in the stride interval of humans during walking of approximately 3-4% vierordt . This random variability has been shown Hausdorff4 ; hausdorff1 ; west1 ; west2 ; griffin  to exhibit long-time correlations, and suggested that the phenomenon of walking is a self-similar, fractal, activity. The existence of fractal time series better suggests that the nonlinear oscillators needed to model locomotion operates in the unstable, that is, in the chaotic regime.

A stochastic version of a CPG was developed by Hausdorff et al. Hausdorff4 ; griffin  to capture the fractal properties of the inter-stride interval time series. This stochastic model was later extended by Ashkenazy et al. ashkenazy  to describe the changing of gait dynamics as we develop from childhood to adulthood. The model is essentially a random walk on a correlated chain, where each node of the chain is a neural center of the kind discussed above, and with a different frequency. This random walk is found to generate a fractal process, with a multifractal width that depends parametrically on the range of the random walker’s step size. Ashkenazy et al. ashkenazy  focused on explaining the changes in the gait time series during maturation, using the stochastic CPG model.

Herein we extend the previous models by assuming that gait dynamics are regulated by a stochastic correlated CPG similar to that of Ashkenazy et al. ashkenazy , coupled to the nonlinear oscillators needed to model locomotion in the unstable, forced and chaotic regimes. We show that two parameters, the average frequency $`f_0`$ and the intensity $`A`$ of the forcing component of the nonlinear oscillator, are sufficient to determine both the fractal and multifractal variability of human gait under normal, stressed and metronomic conditions, using the SCPG model.

## III Human gait analysis

In this section we summarize the main fractal and multifractal characteristics of the stride interval of the human gait data discussed in detail elsewhere scafetta1 . We have analyzed the gait time series of 10 persons in the three different conditions of slow, normal and fast walking. Each time series is approximately one hour long for unconstrained walking, see, for example Fig. 1, for slow, fast and normal walking. Similarly, each time series is approximately 30 minutes long for metronomically constrained walking, see, for example, Fig. 2, for slow, fast and normal walking. Participants in the study had no history of any neuromuscular, respiratory or cardiovascular disorders. They were not taking any medications and had a mean age of 21.7 years (range: 18-29 years); mean height $`1.77\pm 0.08`$ meters and mean weight $`71.8\pm 10.7`$ kg. All subjects provided informed written consent. Subjects walked continuously on level ground around an obstacle-free, long (either 225 or 400 meters), approximately oval path and the stride interval was measured using ultra-thin, force sensitive switches taped inside one shoe. For the metronomic constrained walking, the individuals were told only once, at the beginning of their walk, to synchronize their steps with the metronome. More details regarding the collection of data can be found at Physionet physionet  from where the data were downloaded and in Ref. scafetta1 ; Hausdorff4 .

The fractal and multifractal analysis of the data is done by studying the estimated distribution of the local Hölder exponents using wavelet transforms struzik ; Mallat ; arneodo2 ; Stanley1 . To better understand the meaning of the Hölder exponent $`h`$ we recall that the relation between the Hölder and Hurst exponent $`H`$ 2Mandelbrot  in the continuum limit of a monofratal noise is $`h=H1`$ according to the notation adopted in scafetta1 . Consequently, $`h=0`$ corresponds to pink or 1/f-noise; $`1<h<0.5`$ corresponds to antipersistent noise; $`h=0.5`$ corresponds to uncorrelated Gaussian noise; $`0.5<h<0`$ corresponds to correlated noise; $`h=0.5`$ corresponds to Brownian motion and $`h=1`$ corresponds to black noise. Therefore, the fractal properties of the data can be studied by determining the mean value of the distribution of Hölder exponents and the multifractal properties are given by the width of the distribution itself. However, we stress that a time series of finite length will have a Hölder-exponent distribution with a non-zero width. The existence of such a non-zero width can be a source of confusion between a monofractal time series of finite length and a truly multifractal time series. A multifractal time series can only be distinguished from a monofractal time series of the same length, if the width of its Hölder-exponent distribution is significantly larger than that of a correspondent monofractal time series scafetta1 .

Typical distributions of Hölder exponents, for unconstrained walking of a single individual, are depicted in Fig. 3. Fig. 4 shows the average distributions of the Hölder exponents for the cohort of 10 walkers. Figs. 3 and 4 show that stride interval time series for human gait are characterized by strong persistent fractal properties very close to that of 1/f-noise, $`h0`$. However, normal gait is usually slightly less persistent than both slow and fast gait. The slow gait has the most persistent fluctuations and may present non-stationary properties, $`h>0`$. The slow gait fluctuation may deviate most strongly from person to person. The higher values of the Hölder exponents for both slow and fast gait, relative to normal gait, may be explained as due to a stress condition that increases the persistency and, therefore, the long-time correlation of the fluctuations. Moreover, the regular curves of Fig. 4 show that unconstrained walking is characterized by fractal properties that do not change substantially from one individual to another. Finally, a careful comparison of the widths of the distributions of Hölder exponents for the different gaits with the widths for a corresponding monofractal noise of the same length has proven that the stride interval of human gait is only weakly multifractal scafetta1 . However, the multifractal structure is slightly more prominent for fast and slow gait than for normal gait.

Fig. 5 shows typical distributions of the Hölder exponents for metronome-constrained walking, which is little different from the histograms in Fig. 3. Fig. 6 shows the average distributions of the Hölder exponents for all 10 walkers. The figures clearly indicates that under the constraint of a metronome, the stride interval of human gait become more random and the strong long-time persistence of the 1/f-noise is lost for some individuals. The data present a large variability of the Hölder exponents from persistent to antipersistent fluctuations, that is, the exponent spans the entire range of $`1<h<0`$. However, the metronome constraint usually has a relatively minor effect upon individuals walking normally, the second peak at low Hölder exponents in Fig. 6 being attributable to a single person, who has difficulty with the external cadence. Probably, by walking at a normal speed an individual is more relaxed and he/she walks more naturally. The fast gait appears to be almost uncorrelated noise because the distribution of Hölder exponents is centered close to $`h=0.5`$ characteristic of Gaussian or uncorrelated random noise. Finally, the slow gait presents a large variability from persistent to untipersistent fluctuations.

We notice that some individuals may be unable to walk at a given cadence and their attempts to synchronize the pace results in a continual shifting of the stride interval longer and shorter in the vicinity of an average. For these individuals the phasing is never right and this gives rise to a strong antipersistent signal for all three gait velocities.

In summary, the stride interval of human gait presents a complex behavior that depends on many factors. Walking is a strongly correlated neuronal and biomechanical phenomenon which may be strongly influenced by two different stress mechanisms; (a) a natural stress that increases the correlation of the nervous system that regulates the motion at the changing of the gait regime from a normal relaxed condition, to a consciously forced slower or faster gait regime; (b) a psychological stress due to the constraint of following a fixed external cadence such as a metronome. The metronome breaks the long-time correlation of the natural pace and generates a large fractal variability of the gait regime. In the next section we present a multifractal CPG model able to reproduce these properties.

## IV The SCPG model for human gait

In this section we introduce a model of locomotion that governs the stride interval time series for human gait. As anticipated in the previous sections the model has to simulate a CPG collins1  capable of producing a syncopated correlated output associated with a motocontrol process of the gait cycle. Moreover, the model incorporates two separate and distinct stress mechanisms. One stress mechanism, that has an internal origin, increases the correlation of the time series due to the change in the velocity of the gait from normal to the slower or faster regimes. The second stress mechanism, has an external origin, and decreases the long-time correlation of the time series under the frequency constraint of a metronome. We model this complex phenomenon by assuming that the intensity of the impulses of the firing neural centers regulate only the inner virtual frequency of a forced Van der Pol oscillator holmes . The observed stride interval is assumed to coincide with the actual period of each cycle of the Van der Pol oscillator; a period that depends on the unperturbed inner frequency of the oscillator, the amplitude of the forcing function and the frequency of the forcing function.

Since the frequency of the stepping increases in proportion to the amplitude of the electric stimulation mann , we can assume that the time series of the intensity of the impulses fired by the neural centers is associated to a time series of a virtual frequencies$`\{f_j\}`$. So, in the spirit of the model suggested by Ashkenazy et al. ashkenazy , we assume that the long-time correlated frequency of the SCPG is described by a random walk on a finite-size correlated chain, where each node of the chain is a neural center of the kind discussed above, that fires an impulse with a particular intensity that is associated to a paticular virtual frequency. Ashkenazy et al. ashkenazy  focused on explaining the multifractal changes in the gait time series during maturation from childhood to adulthood, assuming that neural maturation is parametrically associated with the range $`C`$ of the Brownian process that activates the nodes of the finite-size correlated chain of frequencies. Here, we adopt a different approach because we are interested in modeling the gait for human adults operating under different conditions. We assume neural maturation and, therefore, the standard deviation $`C`$ of the random walk process remains constant, whereas the strength of the correlation among the neural centers increases with the change of the velocity of the gait from the normal to the slower or faster regimes. The change of velocity is interpreted as a biological stress. Moreover, contrary to Ashkenazy et al. ashkenazy  we do not add any noise to the output of each node to mimic biological noise. The final output given by the actual frequencies of the gait cycle fluctuates due to the chaotic solutions of the nonlinear oscillators in the SCPG, here that being the forced Van der Pol oscillator. The advantage of using chaos in the model, rather than noise, is that chaos is an intrinsic property of the SCPG dynamics and therefore introduces variability in a controllable way.

We observe that nonlinear oscillators may present chaotic regimes and may be forced by an external frequency holmes , so they may be useful not only to describe the change of phase from the walk, trot, canter and gallop of the quadrupeds, but also the variability of the stride intervals observed in humans. In bipeds it is possible to mimic the movements of the two legs with two nonlinear coupled oscillators. However, because the geometry of the bipeds’ gait, contrary to that for quadrupeds, is unique and the two legs must be shifted by $`\pi `$ rad in phase, we can mimic the biped’s gait with only one nonlinear oscillator. In our model we use a well-known neuronal oscillator model, that is, the forced Van der Pol oscillator collins1 ; holmes  that is defined by the following equation

$$\ddot{x}+\mu \left(x^2p^2\right)\dot{x}+\left(2\pi f_j\right)^2x=A\mathrm{sin}\left(2\pi f_0t\right).$$
(1)
The parameter $`p`$ controls the amplitude of the oscillations, $`\mu `$ controls the degree of non linearity of the oscillator, $`f_j`$ is the inner virtual frequency of the oscillator during the $`jth`$ cycle that is related to the intensity of the the $`jth`$ neural fired impulse, and $`A`$ and $`f_0`$ are respectively the strength and the frequency of the external driver. The frequency of the oscillator would be $`f`$ if $`A=0`$.

We notice that the non linear term, as well as the driver, induce the oscillator to move around a limit cycle. The actual frequency of each cycle may differ from the inner frequency $`f`$. We assume that at the conclusion of each cycle, a new cycle is initiated with a new inner frequency $`f_j`$ produced by the stochastic CPG model while all other parameters are kept constant. However, the simulated stride interval is given by the actual period of each cycle of the Van der Pol oscillator.

We assume that the neural centers of the SCPG may fire impulses with different amplitudes that induce virtual frequencies $`\{f_i\}`$ with finite-size correlations. Here, therefore, we model directly the time series of virtual frequencies. The frequencies $`\{f_i\}`$ are centered around the driver frequency $`f_0`$ according to the relation

$$f_i=f_0+\gamma X_i$$
(2)
where $`\gamma `$ is a constant and $`X_i`$ is a finite-size correlated variable, that is,

$$\frac{<X_iX_{i+r}>}{<X_i^2>}=\mathrm{exp}\left[\frac{r}{r_0}\right].$$
(3)
The parameter $`r_0`$ measures the spatial range of the correlations of the neural network. The chain of frequencies $`f_i=f_0+\gamma X_i`$ is generated by a first-order autoregressive process, also known as a linear Markov process priestley , that is generated by the recursion equation

$$X_i=aX_{i1}+\epsilon _i,$$
(4)
where $`0<a<1`$ is a constant and $`\{\epsilon _i\}`$ is a normalized zero-centered discrete Gaussian process. It is easy to prove priestley  that the autocorrelation function of the chain $`\{X_i\}`$ is given by

$$\rho (r)=\frac{<X_iX_{i+r}>}{<X_i^2>}=a^r.$$
(5)
A direct comparison between Eqs. (3) and (5) gives $`a=\mathrm{exp}[1/r_0]`$, so, we can easily generate a data sequence with the desired finite-size correlation value $`r_0`$. Following Ref. ashkenazy , we assume that a frequency is activated by the position of a random walker given by the function $`g(j)`$ with $`j=1,2,\mathrm{}`$, whose jump sizes follow a Gaussian distribution of width $`C.`$ The width of this distribution, according to the interpretation of Ashkenazy et al. ashkenazy , is associated with the human neural age maturation. This random walk mechanism allows us to obtain from the finite-time, correlated frequency series $`\{f_i\},`$ a new time series of frequencies $`\{f_j\}`$ with $`i=g(j),`$ characterized by long-time correlations. Finally, the sequence of frequencies $`\{f_j\}`$ is used in Eq. (1) recursively.

Normal gait, characterized by the frequency $`f_{0,n}`$, is assumed to occur when the individual is relaxed, consequently the correlations between the neuronal centers are minimum. By implication, whether the gait increases or decreases in velocity, the correlations between the neuronal centers increases. This increase in the stress is modeled by using the short-time correlation parameter $`r_0`$ of the stochastic CPG by assuming

$$r_0=r_{0,n}\left[1+B\left(f_0f_{0,n}\right)^2\right],$$
(6)
where $`r_{0,n}`$ is the short-range correlation among the firing neural centers at the normal frequency gait and $`B`$ is a positive constant that measures the increasing of short-range correlation at the anomalous frequency gait.

We observe that in our experimental data set we measure the stride intervals of the gait. It is true that the frequency of walking may be associated with a long-time correlated neural firing activity that induces virtual pace frequency, nevertheless the walking is also constrained by the biomechanical motocontrol cycle that directly controls movement and produces the pace itself. The SCPG introduced here incorporates both the neural firing activity and the motocontrol constraints.

In summery, the SCPG model is based on the coupling of a stochastic with a hard-wired CPG model and depends on many factors. The most important parameters are the short-correlation size $`r_0`$ of Eq. (3), that measures the correlation between the neuron centers of the stochastic CPG, the intensity $`A`$ of the forcing driving component of the nonlinear oscillator of Eq. (1) and, of course, the average frequency $`f_0`$ of the actual pace that distinguish the slow, normal and fast gait regimes. The other parameters, $`\gamma `$, $`C`$, $`\mu `$ and $`p`$ may be, to a first-order approximation, kept fixed.

While the numerical simulations are left to the following section, we can anticipate an interpretation of the two main parameters $`r_0`$ and $`A`$. In fact, the short-correlation size $`r_0`$ may be interpreted as a parameter that measures the natural correlation between the neural centers and such short-time correlation increases under particular stress, for example, when the velocity of the gait is slower or faster than the normal relaxed situation. The intensity of the forcing driving component $`A`$ may be associated with the voluntary action of trying to follow a particular cadence and is expected to increase under a metronomic constraint.

## V Simulated stride interval gait

In this section we present and comment on our computer simulations of the stride interval of human gait under a variety of conditions. For simplicity, we make use of the following values of the parameters. The frequency of normal gait is fixed at the experimentally determined value of $`f_{0,n}=1/1.1Hz`$, so that the average period of the normal gait is 1.1 second; the frequency of the slow and fast gait are respectively $`f_{0,s}=1/1.45Hz`$ and $`f_{0,f}=1/0.95Hz,`$ with an average period of 1.45 and 0.95 seconds, respectively, that is similar to experimentally realized slow and fast human gaits shown in Fig. 1.

Also the hopping-range parameter is chosen equal to that for adults ashkenazy , that is, $`C=25`$ and kept constant. Moreover, we chose $`r_{0,n}=25`$ such that for $`f_0=f_{0,n}`$ we have $`r_0=25,`$ that coincides with the correspondending value found in Ref. ashkenazy . To generate an artificial sequence with a variability compatible to that of the experimental sequence, we chose $`B=50`$ in Eq. (6) and, in Eq. (2), $`\gamma =0.02`$ , that is a value compatible to the average of the standard deviation of all the data analyzed by us scafetta1 , however, the value of $`\gamma `$ may be smaller and may decrease with an increase in the frequency $`f_0`$ and/or an increase in the intensity of the forcing amplitude $`A`$ of Eq. (1). So, we choose a frequency $`f_0`$, calculate $`r_0`$ via Eq. (6) and the Markovian parameter $`a`$, then we generate a chain of frequencies $`\{f_i\}`$ via Eqs. (2) and (4) and, finally, by using the random walk process to activate a particular frequency of the short-time correlated frequency neural chain, we obtain the time series of the frequencies $`\{f_j\}`$ to use in the time evolution of the Van der Pol oscillator. For simplicity, we keep constant the two parameters of the nonlinear component of the oscillator (1), $`\mu =1`$ and $`p=1`$. The only parameters allowed to change in the model are the mean frequency $`f_0`$ that changes also the value of $`r_0`$ via Eq. (6), and the intensity $`A`$ of the driver of the Van der Pol oscillator (1).

Fig. 7 shows the stride interval time series for slow, normal and fast computer-simulated gaits using SCPG. For the simulation of normal gait we use $`A=1`$ and for both slower and faster gait, we use $`A=2`$. We suppose that the amplitude $`A`$ of the driver of the Van der Pol oscillator (1) should be smaller for normal gait than that for either the slower or faster gaits, because in our interpretation $`A`$ measures the magnitude of the constraint to walk at a particular velocity. The amplitude $`A`$ is smaller for normal gait because normal gait is the most relaxed, spontaneous and, and consequently the most automatic of the three gaits. The figure shows that the SCPG model is able to reproduce a realistic persistence and volatility for the three gaits by simply changing the frequency of the gait itself. In particular note the high volatility of the slow gait that is remarkably similar to that seen in Fig.1.

Fig. 8 shows the stride interval time series for slow, normal and fast metronome-triggered computer-simulated gaits. We use the same frequency series generated by the SCPG used to produce the sequences of Fig. 7. We only change the intensity $`A`$ of the driver of the Van der Pol oscillator (1). We use for the normal gait $`A=4`$ and for both slower and faster gait $`A=8`$. Again we suppose that the intensity $`A`$ of the driver of the Van der Pol oscillator (1) should be smaller than that for both slower and faster gaits, because normal gait is the most relaxed and spontaneous. By comparing Figs. 7 and 8 we note the increase in randomness, the loss of persistency and the reduction in volatility; all effects that are induced in the latter time series by increasing the value of $`A`$ and are found in the phenomenologic data shown in Figs. 1 and 2.

Fig. 9 shows histograms of distributions of the Hölder exponents for the three computer-simulated gaits shown in Fig. 7. The calculation are done in the same way of those used to produce the histograms in Fig. 3 for the experimental data, for details see Ref. scafetta1 . The figure shows that the SCPG model is able to generate artificial stride interval time series with statistical properties similar to the fractal and multifractal behavior of real data. By changing the gait mode from slow to normal, the mean Hölder exponent $`\overline{h}`$ decreases. In the same way by changing the gait mode from normal to fast, the mean Hölder exponent again increases, just as it does for the real data. According to the SCPG model, this increase in the scaling parameter is due to the increase of the inner short-time correlation among the neuronal centers, modeled by Eq. (6). Furthermore, this behavior is due to the biological stress of consciously walking at a speed that is different from the normal spontaneous speed. In addition the multifractality of the gait time series slightly increases for a walking rate different from normal. Here again this effect is observed in the real stride interval data and it is proven by a slight increase in the width of the histograms for fast and in particular slow gait.

Fig. 10 shows the histograms of probability density estimations of the Hölder exponents for the three metronome-triggered computer-simulated gaits shown in Fig. 8. The calculated points show that the SCPG is able to generate artificial stride interval time series that present similar fractal and multifractal behaviors to those of real stride interval data taken under the constraint of a metronome. By increasing the intensity $`A`$ of the driver of the Van der Pol oscillator (1) the randomness of the time series increases and it is possible to obtain a large variety of time series, from those having antipersistent to those with persistent fractal properties. In the SCPG, the parameter $`A`$ measures the constraint of consciousness on the gait, and therefore the value of $`A`$ has to increase if the walker is asked to synchronize his/her pace with the frequency of a metronome. The figure suggests that the SCPG model is able to explain a number of other properties of the metronome-triggered walking. Fig. 6 shows that the usually normal metronome-triggered gait is that with the highest persistent fractal properties. Normal gait is also the most natural under the constraint of the metronome and, therefore, we should expect that the normal gait is the most automatic and the least constrained by human consciousness. This is the reason we have chosen $`A=4`$ for the normal metronome-triggered gait. For both slower and faster metronome-triggered gaits we have chosen $`A=8`$ to indicate a higher conscious stress, that constrains gait at anomalous speeds. Moreover, by comparing Fig. 9 and Fig. 10 and considering that in both simulations we have used the same value of the forcing parameter $`A`$ for both slower and faster gaits, we notice that the largest fractal shift occurs for the slower gait. This increased shift implies that the slower gait is the more sensitive to a voluntary constraint and, so, and so the slower mode has the larger variability. In fact, our human experience and the superposition of the distributions of Hölder exponents for the 10 cohorts in Fig. 6, show a large fractal variability of the slower gait. Finally, Fig. 6 reveals that few persons are characterized by a strong antipersistent pace when asked to follow a metronome. According to the SCPG model , some people are not able to find a natural synchronization and need to continuously adjust and readjust the speed of their pace to match the beat of the metronome. This changing of pace implies a very strong conscious act and, therefore, a very high value of the parameter $`A`$ that would produce a strong antipersistent signal.

## VI Conclusion

We have introduced a new kind of CPG model. One able to mimic the complexity of the stride interval sequences of human gait under the several conditions of slow, normal and fast regimes for both walking freely and keeping the beat of a metronome. The SCPG model is based on the assumption that human locomotion is regulated by both the central nervous system and by a motocontrol system. A network of neurons produces a correlated syncopated output that is correlated according to the level of physiological stress and this network is coupled to the motocontrol process. The combination of systems controls locomotion and the variability of the gait cycle. It is the period of the gait cycle that is measured in the data sets considered herein. Moreover, walking may be conditioned by a voluntary act as well, for example, walking may be consciously forced following the frequency of a metronome. We model the complex system generating the data by assuming that the correlated firing activity of the neural centers regulate only the inner frequency of a forced Van der Pol oscillator. However, it is the forced Van der Pol oscillator that mimics the motocontrol mechanism of the gait cycle. The stride interval is the actual period of each cycle of the forced Van der Pol oscillator. In this way the gait frequency is slightly different from the inner frequency induced by the neural firing activity whose inpulse intensity are able to generate only a potential, but not an actual frequency. The chaotic behavior of such a nonlinear oscillator, like the Van der Pol oscillator, and the possibility to force the frequency of the cycle with an external fixed frequency allows the SCPG model to generate time series that present similar fractal and multifractal properties to that of the human physiological stride interval data in all situations here analyzed. Moreover, by implementing the SCPG with four coupled forced Van der Pol oscillators as in Ref. collins1 , it should be possible to simulate the change of phase between various modes of quadrupeds’ locomotion, that is, walking, trotting, cantering and galloping.

The variety of complex behaviors is regulated by two parameters, the average frequency $`f_0`$ and the amplitude $`A`$ of the driver of the Van der Pol oscillator. The frequency $`f_0`$ regulates the speed and may be associated with a neuronal stress that increases the correlation among the neural centers. The amplitude $`A`$ may be associated with the voluntary action of trying to track a particular frequency and it is expected to increase under a metronome constraint. Finally, Ref. 2ashkenazy  reports that the stride interval time series for elderly subjects and for subjects with Huntington’s diseases are more random than for young healthy subjects. According to the SCPG model, this may be explained by a decrease of the normal short-range correlation among the neural centers that may be associated with a nervous degeneration caused by injury, disease, or aging. This decrease in correlation may be modeled through $`r_{0,n}`$ of Eq. (6). However, the decrease of correlation in the gait of those subjects may also be associated with an increase of the amplitude $`A`$ of the driving force of the Van der Pol oscillator, Eq. (1). In fact, those subjects may also consciously choose to walk more carefully.

——–
Acknowledgment:
N.S. thanks the Army Research Office for support under grant DAAG5598D0002.