text:
# Using Topological Statistics to Detect Determinism in Time Series

## I Introduction

This paper deals with the problem of determining whether the complex behaviour of a single time series may be explained in terms of a deterministic or a stochastic dynamical system. Although the idea was present from the very first days of the nonlinear time series analysis , methods explicitly aiming to detect determinism in time series  have only been published rather recently. The fact that noise can mimic or mask deterministic (e.g. chaotic) behavior in classical measures of chaos (Lyapunov exponents, $`K_2`$ entropy and correlation dimension)  has urged the need for more specific methods to discriminate deterministic from stochastic behaviour.

Measures and vector fields , densities and trajectories , metric and topological , statistical and geometrical, these are, roughly speaking, different terms used to denominate the two broad approaches for investigating the chaotic behaviour of dynamical systems. It is a remarkable fact that most of the aforementioned methods to detect determinism are based in the study of trajectories or vector fields. Considering the noisy character of the inverse problem of time series analysis (from the point of view of nonlinear dynamics), the statistical approach must lie at the heart of any proposed methodology. This becomes relevant as soon as the system underlying a time series is ”excited” beyond the simple bifurcations, so that the geometrical information about the shape of the attractor, or the motion on it, is no longer available. This is also true in the case of stochastic dynamical system, where noise can make impossible to reconstruct the geometrical information . Moreover, trajectories and the geometry of the flow on the attractor are well defined in the case of a deterministic system, in contrast to the case of a stochastic dynamics.

Bearing that in mind, we have developed a novel approach based in the study of the differentiability of the natural measure. The statistical character of our methodology is twofold. The first one concerns the use of the natural (or physical) measure. Differentiability of this invariant measure will be considered in the direction of the evolution, that is, along any typical trajectory. In this way, ”smoothness” of the trajectory and measure are considered in a single step. This is a general result of dissipative systems, regardless of its application to time series analysis. The second one is the use of statistics for analyzing mathematical properties , which allow us to test the differentiability of the measure in a statistical sense. In this way we can extend our methods to signal analysis.
The procedure is based upon a formula that will be derived in the next section, and that will allow us to discriminate between deterministic and stochastic behaviours. This problem was tackled by us in  whereas the basic idea of the methodology was proposed in .

The paper is organized as follows. Section II is devoted to the mathematical background of our approach. At the end of the section we have included a brief summary of the concepts introduced in this section. This is addressed to the less inclined mathematical reader, which can skip most of the mathematical details given in the section. In Sec. II.A we introduce the fundamentals of our approach. The extension to stochastic dynamics is considered in Section II.B. General concepts about measure projections and time series are explained in section II.C, whereas Section II.D is devoted to a brief discussion of the statistics of topological properties used in this work. Section III is devoted to discuss the models to which we apply our method and the numerical procedures followed in each case. We also propose a variant of the surrogation method, partial surrogation, as a way to tune the degree of stochasticity. The results are presented in Section IV. The sensitivity of the time derivative of the measure to stochasticity, the dependence of the results on the embedding dimension and the application of our approach to mixed time series are illustrated on the models discussed in Section III. We also present a simple test on an experimental time series like the Belousov-Zhabotinski chemical reaction.

## II Theoretical Approach

Here we discuss in detail the mathematical background of our method. The natural measure is first defined and the equation of motion that it obeys explicitly given. Then we show how stochasticity triggers a wild behavior of the time derivative of the natural measure. This is the feature upon which our proposal is based. In order to evaluate quantitatively this behavior we borrow the method proposed by Pecora et al . In particular we evaluate (statistically) the continuity of the logarithmic derivative of the natural measure. Those readers more interested in applications may directly go to the Summary at the end of the Section where the most important features of the method are highlighted with the aid of the equations discussed hereafter.

### A The Natural Invariant Measure along the trajectory

Consider a dissipative dynamical system described by $`n`$-first-order differential equations $`\dot{𝐱}=𝐅(𝐱)`$. The corresponding flow $`f^t`$ maps a ”typical” initial condition $`𝐱_0`$ into $`𝐱(t)=f^t(𝐱_0)`$ at time $`t`$. Once transients are over, the motion settles over the attractor $`𝒜`$.

Given such a system, we can define a probability space (measure space) $`(𝒜,,\mu )`$, where $``$ is the $`\sigma `$-algebra generated by the open sets of the invariant set $`𝒜`$, and $`\mu `$ is an invariant measure defined over the sub-sets of $``$ such that $`\mu :(0,1)`$.

Of all the invariant measures which can be defined, only one is relevant from the point of view of the experimentalist or in computer simulations. This is the natural invariant measure, which gives the limiting distribution of almost all starting initial conditions. Using the indicator function $`1_B(𝐲)=1`$, if $`𝐲B_ϵ`$ and $`0`$ otherwise, we can define this measure for a set $`B_ϵ(𝐱)=\{𝐲:d(𝐱,𝐲)ϵ\}`$ as,

$`\mu (B_ϵ(𝐱))=\underset{t\mathrm{}}{lim}{\displaystyle \frac{1}{t}}{\displaystyle _0^t}1_B(f^\tau (𝐲_0))𝑑\tau `$
for almost all $`𝐲_0`$ in the basin of attraction. The so–defined measure is invariant in the sense that it remains constant upon application of the evolution operator $`f^\tau `$,

$`\mu (B_ϵ(𝐱))=\mu (f^\tau (B_ϵ(𝐱)).`$
In the present case what we calculate is $`\mu (B_ϵ(𝐱(t))`$ along the trajectory. This means that the evolution operator changes the point at which the set $`B_ϵ`$ is centered, namely, $`\mu (B_ϵ(f^\tau (𝐱)))`$, leaving unchanged the ball $`B_ϵ`$. This is the so–called Lagrangian evolution (see Fig. 1), that has to be distinguished from the Liouvillian evolution, $`\mu (B_ϵ(𝐱(t))f^\tau (\mu (B_ϵ(𝐱(t))`$, which, contracts the volume along the trajectory, as illustrated in Fig. 1.

Calling $`\mu (𝐱(t))=\mu (B_ϵ(𝐱(𝐭)))`$, the ”material derivative” of the natural measure can be expressed as,

$$\frac{d\mu (𝐱(t))}{dt}=\dot{𝐱}.\mu $$
(1)
or

$$\frac{d\mu (𝐱(t))}{dt}=.(\mu 𝐅)\mu .𝐅$$
(2)
Eq. (2) only involves partial derivatives of $`\mu (𝐱(t))`$ along the trajectory, and since each trajectory in the attractor is contained in the support, $`\mu (𝐱(t))`$ is smooth along any trajectory, and is therefore well-defined. It is known that for axiom-A systems, hyperbolic and with a dense set of periodic orbits, it is possible to decompose the tangent bundle at each $`𝐱𝒜`$ in two linear sub-spaces, stable and unstable. It is also possible to define a measure in these systems, known as a SRB (Sinai-Ruelle-Bowen) which have the property of being smooth along the unstable manifold , due fundamentally to the stretching of the trajectories. On the other side, it is expected to have a wild behavior along the stable direction. In the present case, considering that Eq. (2) only involves derivatives along the trajectory and that in this direction no expansion or compression of the flux occurs (the associated Lyapunov exponent is equal to zero), we expect a smooth behaviour of the measure.

### B Stochastic Dynamics

Consider now a stochastic dissipative dynamical system with additive noise, described by $`n`$-first-order differential equations,

$`\dot{𝐱}=𝐅(𝐱)+\eta 𝐆(t)`$ (3)
where $`\eta >0`$ is a small number (noise intensity) and $`𝐆(t)`$ is a vector of independently and identically distributed random Gaussian variables, of zero mean and correlations $`<G_i(t)G_j(t^{})>=\delta _{ij}\delta (tt^{})`$. A physical system will normally have a small level $`\eta `$ of random noise, so that it can be considered a stochastic process rather than a deterministic one. In a computer study, roundoff errors should play the role of the random noise. For suitable noise and $`\eta `$, the stochastic time evolution (3) has a unique stationary measure $`\mu `$ . This is the natural (or physical) invariant measure defined above.

Inserting Eq. (3) in (1),

$$\frac{d\mu (𝐱(t))}{dt}=(𝐅+\eta 𝐆(t)).\mu $$
(4)
Whenever the vector field $`𝐅(𝐱)`$ can be expressed as $`𝐅(𝐱)=\varphi (𝐱)+𝐟(𝐱)`$, with $`𝐟(𝐱)`$ being orthogonal to the gradient term and having no divergence, the measure is given by :

$$\mu (𝐱)=N\mathrm{exp}(\frac{\varphi (𝐱)}{\eta ^2})$$
(5)
Then, introducing (5) in (4) we arrive at:

$$\frac{d(\mathrm{ln}\mu (𝐱(t)))}{dt}=\frac{1}{\eta }(\frac{1}{\eta }|𝐅(𝐱)|^2+𝐆(t).𝐅(𝐱))$$
(6)
In section II.A we discussed the smoothness of the measure along the trajectory for a deterministic system, at least in the case there exists a SRB measure. In the present case, that is, a stochastic process, it is possible to define a unique stationary measure that will tend to the SRB for $`\eta 0`$. This is the Kolmogorov Measure . At each time step of evolution, we add some noise (of amplitude $`\eta `$). If we repeat these operations (evolving by time evolution and putting some noise) again and again then we will get an invariant measure which is smooth along the unstable direction. This is because the deterministic part of the time evolution will improve the continuity of the density in the unstable direction by stretching, and roughening it in the other directions due to contraction. Thus, in the zero-noise limit $`\eta 0`$ we will get a measure $`\mu `$ that satisfies SRB conditions.

Eq. (6) provide an alternative tool to investigate the findings of , according to which smoothness in phase space implies determinism in time series. For weak noise levels, the first term of the right hand side of (6) is dominant over the second term. In this case, smoothness in phase space implies ”continuity” in the left hand side of (6), or differentiability of the measure along the trajectory. On the other hand, in the case of strong noise levels, the second term is the dominant one and a wild behavior in the measure must be expected. This is so because the vector $`𝐆(t)`$ is uncorrelated with the actual position in the phase space. Although Eq. (6) was derived for a vector field obeying some restrictions (see above), our numerical results indicate that it can be applied to more general systems.

### C Time Series and Measure Projection

In case that dynamical invariants are to be estimated from the knowledge of a single observable, we can rely upon Takens’ Delay Coordinate Map Theorem . According to it, the structure of the attractor (including differential information) can be preserved by using delay coordinate maps , as long as the embedding dimension $`m`$ is large enough to fully unfold the attractor structure. Furthermore, the Fractal Whitney Embedding Prevalence Theorem  tells us that almost every smooth map will be an embedding (one-to-one and differential structure preserving) provided that $`m>2D`$, where $`D`$ is the box-counting dimension of the attractor. However, as shown in Refs. , when dealing with functions of the measure, embedding dimensions greater than the correlation dimension are enough.

From the statistical point of view, we must ask how is the embedding process in the measure space $`(A,)`$ associated to the attractor $`A`$. That is, if $`\mu `$ is a measure in $`\mathrm{}^n`$ and $`\varphi :\mathrm{}^n\mathrm{}^m`$ is a function, we will have a projected or induced measure $`\mu _\pi =\varphi (\mu )`$ over a subset $`S\mathrm{}^m`$, where $`\varphi (\mu )(S)=\mu (\varphi ^1(S))`$. This is important because we are actually transforming the measure in the embedding process. Most of the existing literature is related with projections and reconstructions of the invariants sets and the dimensions of the corresponding probability measures . In our case, if $`\mu `$ is the physical measure describing the original system ($`\mu `$ is carried by an attractor in phase space) then the points in the reconstructed space are equidistributed with respect to the projected measure $`\mu _\pi `$  (except for particular cases avoided explicitly by the embedding theorems ). This is a sufficient condition which allows us to extract information about the original system working in its corresponding projection.

If we have a time series, we can construct an $`m`$-dimensional vector:

$`\overline{x_i}=(x_i,x_{i+\tau },\mathrm{},x_{i+(m1).\tau })`$
where $`\tau `$ is the time delay, chosen by one of the standard methods . Then, from a single set of observations, multivariate vectors in $`m`$-dimensional space are used to trace out the orbits of the system. Fig. 2 shows the concept of the reconstruction method. Explicitly shown is the reconstruction of the measure carried by the attractor $`\mu _\pi `$.

### D Statistics of Topological Properties

Taken’s theorem , and its sequels, give a rigorous justification for state space reconstruction. The essence of the mathematical proof is that the trajectory formed from the time series is diffeomorphically related to the actual phase-space trajectory of the dynamical system. In order to test the mathematical properties embodied in the diffeomorfism, that is, continuity, differentiability, inverse differentiability and injectivity, Pecora et al.  have developed a set of statistics intended to test quantitatively these properties. Their algorithms are of general use and can in particular be applied to test topological properties in any pair of sets of points. We will take advantage of this procedure by implementing the aforementioned algorithms to test numerically the continuity of the logarithmic derivative of the measure along the trajectory. Basically, the method is intended to evaluate, in terms of probability or confidence levels, whether two data sets are related by a mapping having the continuity property: A function $`f`$ is said to be continuous at a point $`𝐱_0`$ if $`ϵ>0,\delta >0`$ such that $`𝐱𝐱_0<\delta f(𝐱))f(𝐱_0)<ϵ`$. The results are tested against the null–hypothesis, specifically, the case in which no functional relation between points along the trajectory and the measure exists. This is done by means of the statistics proposed by Pecora et al

$`\mathrm{\Theta }_{C^0}(ϵ)={\displaystyle \frac{1}{n_p}}{\displaystyle \underset{j=1}{\overset{n_p}{}}}\mathrm{\Theta }_{C^0}(ϵ,j)`$ (7)
and

$`\mathrm{\Theta }_{C^0}(ϵ,j)=1{\displaystyle \frac{p_j}{p_{max}}}`$ (8)
where $`p_j`$ is the probability that all of the points in the $`\delta `$-set, around the point $`𝐱_j𝐱(t)`$, fall in the $`ϵ`$-set around $`\frac{d\mathrm{ln}\mu (𝐱_j)}{dt}`$. The likelihood that this will happen must be relative to the most likely event under the null hypothesis, $`p_{max}`$ (see reference ). When $`\mathrm{\Theta }_{C^0}(ϵ,j)1`$ we can confidently reject the null hypothesis, and assume that there exists a continuous function. As in the work of Pecora et al  the $`ϵ`$ scale is relative to the standard deviation of the density time series, and thus, $`ϵ[0,1]`$. Plots of $`\mathrm{\Theta }_{C^0}(ϵ)`$ versus $`ϵ`$ can be used to quantify the degree of statistical continuity of a given function. In order to characterize the continuity statistics by means of a single parameter we have also calculated,

$$\theta =_0^1\mathrm{\Theta }_{C^0}(ϵ)𝑑ϵ$$
(9)
The limiting values of $`\theta `$, namely, 0 and 1, correspond to a strongly discontinuous and a fully continuous function, respectively.

### E Summary

A naive test to quantify noise in signals is to check how smooth they are. As long as more noise contaminates the signal, more discontinuous it becomes. This is the case for example of additive noise, e.g. noise added to the signal. However, this is by no means a general rule. For instance, intrinsic noise, that is, noise added in the equation of motion, is not expected to affect the smoothness of the signal. This can be clearly seen in the case of surrogate time series: two time series (the original series and its surrogate) having the same correlation structure, may have very different underlying dynamics, i.e. one deterministic and the other stochastic. We can overcome the above drawback by using the trajectory of the system, instead of a single variable. Smoothness or continuity of the trajectory in phase space has been used before  in this context. What we propose here is to use the distribution of points on the trajectory (or the natural measure) as a way to evaluate the the degree of noise in the system. Our proposal is based upon Eq. (6). This equation tells us that the amount of noise present in the system is directly related to the differentiability of the natural measure, evaluated along a typical trajectory. More noisy is the system, less differentiable it becomes. It remains the question of how to apply a fundamental analysis concept, like differentiability, to a magnitude defined numerically. The ”topological statistics” tools developed by Pecora et al  came to our help. It is devised to evaluate in a statistical fashion topological properties of functions. Eq. (8) is the statistics used to test, against the null hypothesis, the degree of continuity of our function. As long as this statistic approaches unity we are confident that we have a continuous function. We must note that we have evaluated the continuity of the numerical derivative of the measure, which is equivalent to test differentiability, a computationally easier procedure. Continuity is directly related to the resolution with which we are looking at the function, that is $`ϵ`$, so the statistics is actually $`ϵ`$–dependent $`\mathrm{\Theta }_{C^0}(ϵ)`$. In order to provide us with an overall continuity test, we summed up over the whole range of $`ϵ`$ obtaining a single parameter $`\theta `$ (see Eq. (9)) which we use hereafter to characterize continuity.

## III Models and Numerical Procedures

Here we discuss the various stochastic dynamical systems on which we have investigated the efficiency of our approach. The numerical procedures followed to reconstruct the space from time series related to some of the coordinates of those models and to evaluate the statistical differentiability of the natural measure are also described. We also discuss a variant of the surrogation process which consists of producing partially surrogate series. The method is a way to vary the degree of stochasticity of the time series.

### A The stochastic Van der Pol Oscillator

The simplest case in which we can apply our ideas is in a nonlinear oscillator. We have used the Van der Pol (VdP) oscillator  with an additive stochastic term.

$`\dot{x}`$ $`=`$ $`y+\eta G_1`$ (10)
$`\dot{y}`$ $`=`$ $`(x^21)yx+\eta G_2`$ (11)
The parameter $`\eta `$ represents the noise level, and $`G_i(t)`$ are uncorrelated Gaussian noises, such that $`G_i(t)`$ Normal(0,$`\sigma `$), zero mean and standard deviation $`\sigma `$. Without loss of generality we take $`\eta =1`$, and tune the degree of noise only by the standard deviation $`\sigma `$.

Numerical integration was carried out by means of an Euler algorithm. 20000 data points have been generated. Using the $`x`$-coordinate as our “experimental” time series, we have made a reconstruction with an embedding dimension in the range 2–20, and a $`\tau `$ lag of 10 (in sampling units). Every reconstruction has been rescaled to the unit hyper-square.

Density time series have been obtained for each reconstruction. Starting with the first point in the reconstructed phase space, we follow the trajectory recording the density of points around each of the trajectory’s point. In order to estimate this density, we have used the Epanechnikov kernel , which, roughly speaking,“weighs” the points according to its distance to the center. This is preferable to the Gaussian kernel, which is of infinite support (and therefore has a lower computational efficiency) and the “square” kernel which gives equal status to the different points in the ball around the reference point. A parameter that has to be chosen carefully is the ball radius used to estimate the density. If the radius is too small, the low density regions will be practically depopulated and the measure will be underestimated. On the other hand, if the radius is too big, the estimation will capture points which are not really part of the neighborhood of the reference point. Of course, the ball radius is a function of the data points being used in the reconstruction process and must be chosen according to this fact. We have used a radius of the ball of 5% of the total attractor extent. In evaluating the continuity statistics, we average $`\mathrm{\Theta }_{C^0}(ϵ,j)`$ over $`n_p`$ points (see Eq. (6)) randomly distributed in the trajectory, typically $`10\%`$ of the total record.

Fig. 3 shows a typical density time series from the $`x`$-coordinate of the Van der Pol oscillator. Albeit qualitative, the smooth behavior of this density along the trajectory is readily noted.

### B The stochastic Lorenz System

The VdP oscillator will allow us to get a closer look at the procedure we are implementing. However, a nonlinear oscillator is a somewhat trivial example and we want to apply the method to more complicated cases. In fact, the ultimate objective of the methodology is to discriminate random behaviour from a deterministic one, and this is specially important in the case of chaotic behaviour.

The Lorenz system  with an additive stochastic term is an adequate choice. The related system of differential equations can be written as:

$`\dot{x}`$ $`=`$ $`sx+sy+\eta G_1`$ (13)
$`\dot{y}`$ $`=`$ $`y+rxxz+\eta G_2`$ (14)
$`\dot{z}`$ $`=`$ $`bz+xy+\eta G_3`$ (15)
The parameters used in the calculations are, $`s`$ = 10.0, $`r`$ = 28.0 y $`b`$ = 2.66, which give chaotic behaviour in the case $`\eta =0`$. $`G_i(t)`$ and $`\eta `$ were defined above. When not specified the results for the Lorenz system discussed in the following Section were obtained for an embedding dimension of 3, which is greater than the correlation dimension of the Lorenz system.

Numerical integration of the Lorenz system was carried out by means of the Euler method. The time integration step was 0.01. Time series with 16384 data points and their respective surrogates were subsequently generated. The reconstruction was performed by the usual time–delay method , with a time delay given by the first zero of the autocorrelation estimate (10 in units of the integration step) on an embedded phase space of dimension 3. The natural measure $`\mu (𝐱(t))`$ along the trajectory was calculated by means of the Epanechnikov kernel density estimator  with an sphere of radius 5% of the attractor extent. As in the VdP oscillator the continuity statistics was evaluated including up to 10% of the points in a given record.

The sensitivity of the time derivative of the measure to stochasticity is illustrated in Fig. 4. This figure shows that whereas the surrogate of the $`x`$–coordinate time series is as ”smooth” as the original series, the time derivative of the logarithm of the measure is much more spiked in the surrogate than in the original series.

### C The Mackey–Glass model

In order to investigate the effects of the embedding dimension we have also considered the high dimensional system introduced in . The dynamical system is described by means of the following delay–differential equation,

$$\dot{x}=\frac{ax(t\delta }{1+x(t\delta )^c}bx(t)$$
(16)
This equation has been proposed to model nonlinear feedback control in physiology. We use the set of parameters that gives an attractor dimension of $`7.5`$ , namely, $`a=0.2`$, $`b=0.1`$, $`c=10`$, and $`\delta =100`$. The resulting time series were analyzed with a time delay given by the first zero of the autocorrelation estimate, and the measure was evaluated on spheres of radius 10% of the attractor extent and the continuity statistics with 10% of the points in the total record.

### D Partial surrogation

The surrogation process is a well established method in the context of nonlinear time series analysis . In typical applications a single time series is available. From this time series, an ensemble of the so-called surrogate series are generated that mimic certain properties of the original. For example, by simply scrambling the temporal order of the points in the original, one obtain surrogate time series which preserve the mean, variance, etc. One of the most popular methods of producing surrogate time series consists of shuffling the phases in the Fourier transform of the original data set . In this way, each value of the Fourier transform of the original data is multiplied by a random phase $`\mathrm{exp}(i\varphi )`$, with $`\varphi `$ random $`[0,2\pi ]`$ <sup>*</sup><sup>*</sup>*In order to get a real time series in the antitransformation we multiply symmetrically with respect to the center of the transform. The procedure generates a new time series with the autocorrelation structure of the original. Here we propose to introduce an additional factor, $`\mathrm{exp}(i\varphi \alpha )`$, with $`\alpha [0,1]`$. This factor allows us to control the degree of stochasticity by tuning the parameter $`\alpha `$.

## IV Results

### A Sensitivity of the time derivative of the measure to stochasticity

In order to test the efficiency of the approach proposed here we have first evaluated the continuity statistics either on a coordinate or on the time derivative of the natural measure time series of the stochastic Van der Pol oscillator. Fig. 5 shows the continuity statistics for the $`x`$–coordinate and for the time derivative of the reconstructed natural measure with and without noise. There is almost no modification in the statistics of the coordinate upon noise addition. However, the statistics of the density time series reflects very clearly the presence of the stochastic term. This illustrates the efficiency and novelty of our approach and supports its application to more complex cases. We must remark that we are using the commonly named ”dynamical” noise, that is, a stochastic term added in the dynamical equations, instead of the ”measurement” noise, which is added after the ”clean” integration step. Preliminary results shows that an efficient method to discriminate both types of noise, can be achieved by using the continuity statistics over the density and the coordinates, but this deserve further research, and eventually will be published elsewhere.

The results for the continuity statistics of the time derivative of the reconstructed measure from the $`x`$–coordinate of the Lorenz system are illustrated in Fig. 6. The results of Fig. 6a show that the time derivative of the measure in the original series is ”more continuous” (in a statistical sense) than its surrogate. Partial surrogation ($`10\%`$) decreases the degree of continuity of the time derivative of the measure in an extent lower than total surrogation, as expected. On the other hand, the results show that the continuity of the totally surrogated series shows almost no dependence on whether it has been derived from the original series or from a partially surrogated series (10 % surrogated).

The results for the stochastic Lorenz system reported in Fig. 6b clearly show that the stochastic terms significantly decrease the statistical continuity of the time derivative of the measure. Surrogation of the stochastic series produces a further decrease of continuity, indicating that the series still has some degree of determinism. The degree of stochasticity of a time series can be quantified by calculating the integral of the continuity statistics as defined in Eq. (8). Fig. 7 a) shows how steeply $`\theta `$ decreases with the percentage of surrogation. Similarly, $`\theta `$ decreases with the standard deviation of the Gaussian noise in the stochastic Lorenz system (Fig. 7b), as expected. Thus, the magnitude $`\theta `$ can be used to evaluate the relative stochasticities of a set of experimental time series.

### B Dependence on the Embedding Dimension

A point of crucial relevance is how the above results change with the embedding dimension $`m`$. We have investigated this question on the Lorenz system and on the high–dimensional system discussed in III.D . The results for the Lorenz system depicted in Fig. 8a show that $`\theta `$ decreases with the embedding dimension. This is a consequence of working with a fixed sphere radius for all $`m`$ and of the numerical noise that should increase with $`m`$. The decrease of $`\theta `$ is stronger in the surrogate series, although it is likely that the difference between the two should decrease for large enough $`m`$. In any case, the difference in $`\theta `$ between the original and the surrogate series changes only from 0.35 to 0.49 when $`m`$ is varied in the range 3–10. The behavior of $`\theta `$ in the high–dimensional system is far more intricate (see Fig. 8b). For $`m`$ well below the attractor dimension the measure for the surrogate series seems to be more continuous than that for the original series. The reason for this rather odd behavior has to be found in the heavy crossing of trajectories that occur at $`m`$ far below the attractor dimension . In those cases, surrogation seems to have a smoothing effect. Instead, for $`m>7`$ the behavior is similar to that of the Lorenz system, although the difference between the original and the surrogated records is substantially smaller. This point further study that is actually in progress.

Fig. 9 illustrates the dependence of the continuity statistics of the measure on the embedding dimension $`m`$ for the stochastic Van der Pol oscillator. In this case the simplicity of the attractor gives an almost null dependence on $`m`$ for no noise or very low noise levels. When the noise is increased the behavior is in line with that found in the Lorenz system, namely, a decrease of $`\theta `$ as $`m`$ increases. This dependence on $`m`$ is more noticeable the greater the noise level.

The study of the dependence of the continuity statistics of the measure derivative on the embedding dimension involves numerical difficulties that deserve some comments. In the case of low enough embedding dimension and moderate number of data points, the estimation can be achieved confidently for almost all the points in the trajectory. Of course, ”moderate” and ”long enough” are terms which depend on the attractor dimension. We think that the same criteria followed to get for example a reliable estimation of the correlation dimension , apply in this case. In our approach we have introduced noise explicitly, adding another factor to be considered in the estimation step. As it is well known, as ”extra” dimensions are included in the embedding process, noise populates them more or less uniformly. This is specially problematic in the case of ”simple” systems, like the VdP oscillator, as the trajectory in a stochastic oscillator may ”wander” far from the zero-noise limit cycle. In these excursions, the measure swept by the trajectory is unavoidably constant, because no other points are in the neiborhood of the evolution, except those which are time correlated. In such a case, a constant measure results, and thus, a high value of the continuity statistics. This effect is more noticeable as the embedding dimension is increased, because ”more” space is available (see Fig. 10).

### C Application to Mixed Time Series

A distinctive feature of our method is the possibility of using it in different ranges of a given time series. In this way we can examine short records and evaluate its stochasticity. Bearing this in mind we have devised the following example: Suppose we have a time series which is half deterministic and half stochastic. Could our method discriminate both behaviors in the same time series?. In order to answer this question, we have generated a single time series (16384 points) with the first half coming from the $`x`$-coordinate of the deterministic Lorenz system, and the second half coming from its surrogate (100 % randomization) time series. We have applied the continuity statistics over four regions in the time derivative of the density record (two randomly selected in the first half and two in the second half). Fig. 11 shows the results. It is then clear that the statistics utilized here can discriminate stochastic from deterministic behaviour. Fig. 11 also shows the statistic for the whole time series (same number of reference points randomly selected along the time series). The results are midway between those for the stochastic and deterministic ranges.

As another example of the use of our methodology we have investigated the effect of a noise burst in the Van der Pol oscillator. The dynamic equations are rewritten as,

$`\dot{x}`$ $`=`$ $`y+A(t)\eta G_1`$ (17)
$`\dot{y}`$ $`=`$ $`(x^21)yx+A(t)\eta G_2`$ (18)
where $`A(t)=1`$ if $`t_1<t<t_2`$ and $`0`$ otherwise. The interval $`[t_1,t_2]`$ is a small interval where we ”turn on” the stochastic term. The idea behind this system is to test the capabilities of the method to detect the noise introduced. Using 20000 data points, we have used 500 consecutive points with the stochastic term added. In Fig. 12 we show a typical time series (and the measure time derivative) where noise has been turned on in the time interval 10000–10500, with a strength of $`\sigma =0.05`$. Fig. 13 shows Pecora et al.  statistics applied to five regions of the whole series, one of them being the stochastic region. Again, our method clearly discriminates noisy and clean regions.

### D Application to Experimental Time Series

In order to test our method in actual experimental time series, we have used the data set from the Belusov-Zhabotinski (BZ) chemical reaction . As shown by those authors the apparently random behaviour of the amplitude of the concentration of bromide ions can in fact be explained by deterministic laws. We will use here our approach to confirm the above finding.

Using the bromide concentration time series, as in the work of Roux et al.  we repeat the procedure explained above. We have used a $`\tau `$ = 30 (in sampling units) for the reconstruction and an embedding dimension of 4. Fig. 14 shows the continuity statistics for both the time derivative of the reconstructed measure using the BZ time series and for its 100 % surrogate. The large difference between the two and the rather high value of $`\mathrm{\Theta }_{C^0}`$ confirm our proposal in the sense that one can explain the behaviour of this record as the output of a deterministic dynamical system.

## V Concluding Remarks

In brief, we have proposed a method to identify determinism in time series which exploits the continuity of the logarithmic time derivative of the natural measure along the trajectory, that is, its differentiability. The method is based upon a formula which explicitly shows the sensitivity of the measure to stochasticity. In the present work we have adapted the statistical method of Pecora et al.  to investigate the continuity of the time derivative of the measure. Results of partially surrogated series and series derived from two stochastic dynamical systems and a high–dimensional system, clearly illustrate the suitability of the present method to the problem at hand. As we have shown, the method is very easily applied to any kind of time series, from the simplest one, as is the case of limit cycle oscillators, to the high–dimensional cases. Given a time series, and once done the standard reconstruction process, obtaining the density (or its time derivative) is a straightforward procedure. Then, one can use Pecora et al. approach over the whole time series or in selected pieces of it. A warning against a blind application of the method is in order, particularly in what concerns the evaluation of the continuity of the measure time derivative. A careful inspection of the density time series must be done before any further operation is performed.

The dependence of the continuity statistics on the embedding dimension in low and high–dimensional systems, indicate that applications to real (experimental) time series would eventually require a thorough investigation of this point in each particular case, as is common in time series analysis. In any case, the fact that the method works reasonably well on short time series, supports its usefulness for the analysis of experimental series. We have shown this in a simple case as it is the Belousov-Zhabotinski reaction, confirming previous findings.

In a broader sense, the application of continuity statistics over the density time series is a new aspect of the possibilities offered by the Lagrangian Measures . As we have shown previously, the use of more traditional tools on this density, such as Fourier Transforms or histograms, may help in extracting new information on the underlying dynamical system.

###### Acknowledgements.

Thanks are due to E. Hernández, M. SanMiguel and R. Toral for many useful comments and suggestions. This work was supported by grants of the spanish CICYT (grant no. PB96–0085), the European TMR Network-Fractals c.n. FMRXCT980183, and of the Universidad Nacional de Quilmes. G. Ortega is member of CONICET Argentina and also thanks the ”Generalitat Valenciana” for support.


format: text

text:
# Does the quark-gluon plasma contain stable hadronic bubbles?

## I Introduction

In a previous paper  we calculated the free energy of a bubble of hadrons embedded in an extended quark-gluon plasma (QGP), and of a droplet of QGP embedded in an extended hadron phase, for parameters in the vicinity of the cosmological quark-hadron transition, i.e. the baryon chemical potential were set to zero. In Ref.  we found, as Mardor and Svetitsky in , that the free energy of a hadron bubble of radius $`R`$ embedded in QGP possessed a minimum at radii of a few fm, even above the phase transition temperature. Thus, hadronization from QGP is strongly enhanced compared to the usual nucleation scenario, where an energy barrier has to be overcome before bubbles of hadrons can form and grow.

An important ingredient in these calculations is the density of states of the relevant particles. In Ref.  we advised a modification of the usual expressions arising from the multiple reflection expansion (MRE), and we saw that this modification (the MMRE) yielded more accurately the free energy than the MRE did when compared to a direct (but numerically heavier) sum-over-discrete-states calculation. Our modification of the multiple reflection expansion expressions is physically well motivated, in that it consists of a truncation of the density of states in such a way that we avoid the effects of a negative density of states, present in the usual MRE expressions.

In this paper, we generalize the calculations to finite chemical potential, and we shall see that also in the case of non-zero chemical potential, the MMRE produces accurately the thermodynamic potential. We show that, as in the case of zero chemical potential, the thermodynamic potential of a hadron bubble embedded in QGP has a minimum at a radius of $`R1`$ fm, meaning that it is thermodynamically favorable for a QGP in equilibrium to spontaneously create bubbles of hadrons of about this size. However, the minimum in the thermodynamic potential found for non-zero chemical potential is less pronounced than in the zero chemical potential case considered in Ref. . The more general calculations presented in this paper should be relevant also to neutron star phenomenology as well as ultrarelativistic heavy ion collisions.

In the following section we give an overview of the basic theory leading to the results of section III. Finally in Sec. IV, we summarize the conclusions.

## II Theoretical framework

### A The particle model

#### 1 The quark-gluon plasma

The model of the plasma phase is a close derivative of the MIT Bag model . It consists essentially of non-interacting quarks and gluons (also self-interaction of the gluons is neglected), bounded by MIT Bag-like boundary conditions. There are three quark flavors in our model: $`u`$ and $`d`$, which we consider massless, and $`s`$, to which we assign a mass of $`m_s=150`$ MeV. The other quarks are too heavy to be relevant in this analysis. We choose a critical temperature (the bulk phase equilibrium temperature at zero chemical potential) of $`T=150`$ MeV, which gives a Bag constant of $`B=(221`$ MeV$`)^4`$. The plasma may be inside a sphere, as in the usual MIT Bag (we shall refer to this configuration as a plasma droplet), or outside, i.e. in a “vacuum bubble” configuration. For a detailed review of this model, the reader may consult . We are mostly interested in the vacuum bubble situation, where the vacuum bubble is filled with a thermal mixture of hadrons as described later.

For the density of states of the quarks and gluons we use a modification of the multiple reflection expansion (MRE) , advised in . This modification affects only the gluon density of states, which is then

$$\rho _g(k,R)=\{\begin{array}{cc}0,\hfill & 0kR<0.832\hfill \\ \frac{V_{QGP}k^2}{2\pi ^2}\frac{4R}{3\pi },\hfill & kR0.832\hfill \end{array}$$
(1)
for gluons inside a sphere (i.e. the MMRE($`R`$)), and

$$\rho _g(k,R)=\{\begin{array}{cc}0,\hfill & 0kR<0.458\hfill \\ \frac{V_Hk^2}{2\pi ^2}+\frac{4R}{3\pi },\hfill & kR0.458\hfill \end{array}$$
(2)
for the difference between the gluon density of states in a volume $`V_{\mathrm{}}`$ and in volume $`V_{\mathrm{}}V_H`$, (i.e. the MMRE($`R`$)).

In  we saw that this modified multiple reflection expansion (MMRE) description of the density of states made the free energy agree nicely with more direct calculations, i.e. summing over discrete energy levels. Before using the MMRE in the more general case under consideration in this paper, we performed similar checks of the MMRE against direct sum-over-discrete-states calculations. Also in this case, the MMRE proved to describe the density of states adequately and significantly better than the MRE, cf. Fig. 1. For a complete description of the MMRE density of states, we refer the reader to .

#### 2 The hadron phase

We include all hadrons with masses below $`1.2`$ GeV, taking only the volume part of the density of states into account. This is justified as the hadrons represent far fewer degrees of freedom than the QGP. Specifically, the density of states of a hadron species occupying a volume $`V_H=\frac{4\pi }{3}R^3`$ is

$$\rho _H(k,R)=\frac{V_Hk^2}{2\pi ^2}.$$
(3)
### B Thermodynamics

The thermodynamic potential is

$$\mathrm{\Omega }(T,V,\{\mu _i\})=T\mathrm{ln}(Z(T,V,\{\mu _i\})),$$
(4)
where $`Z(T,V,\{\mu _i\})`$ is the partition function

$$Z(T,V,\{\mu _i\})=Tr\left\{e^{(_i\mu _i\mathrm{\Lambda }_i)/T}\right\}.$$
(5)
Here, $``$ is the Hamiltonian, the $`\mathrm{\Lambda }_i`$’s are the symmetries of the Hamiltonian (the conserved quantities), and $`\{\mu _i\}`$ denotes the corresponding collection of chemical potentials. $`Tr\{..\}`$ means the sum over all diagonal elements of energy- and number-eigenstates. These are simultaneous eigenstates since $`[,\mathrm{\Lambda }_i]=0`$. The volume dependence enters via the energy levels, i.e. in the trace operation.

In thermodynamic equilibrium, the configuration realized in Nature is the one which minimizes the thermodynamic potential, leading to the Gibbs conditions for phase equilibrium: $`T_1=T_2`$ (thermal equilibrium), $`\frac{\mathrm{\Omega }_1}{V_1}=\frac{\mathrm{\Omega }_2}{V_2}`$ (mechanical equilibrium), and $`\mu _i^{(1)}=\mu _i^{(2)}`$ (chemical equilibrium), where index $`1`$ and $`2`$ refer to the two phases.

If there are no interactions between the particles, i.e. the energy of a particular state is the sum of the energies of the individual particles, then we can write the thermodynamic potential for a particular particle species enclosed in a spherical volume of radius $`R`$, as

$$\mathrm{\Omega }(R,T,\mu )=gT_0^{\mathrm{}}𝑑k\rho (k,R)\mathrm{ln}(1\pm e^{(\sqrt{k^2+m^2}\mu )/T}),$$
(6)
(upper sign for fermions, lower sign for bosons) where we further assume that the energy levels are sufficiently closely spaced that we may use a smoothed density of states, $`\rho (k,R)`$, (e.g. the MMRE) normalized such that $`_0^{\mathrm{}}𝑑k\rho (k,R)`$ is the total number of states in the volume $`V=\frac{4\pi }{3}R^3`$. The chemical potential entering in (6), is the combined chemical potential of that particle species

$$\mu =\underset{i}{}\lambda _i\mu _i,$$
(7)
$`\lambda _i`$ being the quantum expectation value of the symmetry operator $`\mathrm{\Lambda }_i`$ of the particle species in question (e.g. $`\lambda _{\mathrm{b}aryon}=1/3`$ for a quark). Finally, $`g`$ is an appropriate degeneracy factor.

At the energies relevant in this analysis, the quarks and gluons can be considered point-like. This is not the case for the hadrons. We must take explicit account of the fact that the hadrons occupy a certain volume. A number of such excluded volume corrections have been discussed in the literature, but we shall not go into an analysis of the different suggestions. We choose the one proposed in , since it is very easy to implement and grasps much of the essential physics involved. According to , the “true” pressure of the hadron gas ($`p_H`$) is obtained from the ideal gas pressure ($`p_{H,id}`$) according to

$$p_H=\frac{p_{H,id}}{1+ϵ_{H,id}/(4B)},$$
(8)
where

$$ϵ_{H,id}=\underset{i}{}ϵ_{i,id}$$
(9)
is the total energy density of the ideal hadron gas,

$$p_{H,id}=\underset{i}{}p_{i,id}$$
(10)
is the total hadronic ideal gas pressure, ($`i\mathrm{h}adronspecies`$), and $`B`$ is the Bag constant. We note that some sort of excluded volume correction is essential, since in a point-like hadron gas, no transition to QGP will occur at low temperatures, even at arbitrarily high density.

There is one more question about the hadrons we need to address before we proceed. This is the phenomenon of Bose-Einstein condensation, regarding the bosonic part of the hadron gas. A bosonic gas always has $`\mu _im_i`$, and when $`\mu _i=m_i`$ particles will start to “condense” into the lowest, zero-momentum state, meaning that there will be a macroscopic number of particles in this state. When the parameters are such that we would have $`\mu _im_i`$, the hadron species $`i`$ is removed from the above summations (9) and (10), since particles in a zero-momentum state do not contribute any pressure.

## III Results

We shall consider both the case of equal chemical potentials, $`\mu _u=\mu _d=\mu _s\mu _q`$, and the case where the $`u`$ and $`d`$ quarks have a different chemical potential from that of the $`s`$ quark, $`\mu _u=\mu _d\mu _s`$. The first case is relevant when the time scales are such that the weak interactions maintain equilibrium between all three flavors of quarks (e.g. quark matter in neutron stars, neglecting the electron chemical potential as a first approximation), whereas the case of separate $`s`$ quark chemical potential is relevant to ultrarelativistic heavy ion collisions, the time scales here being such that the net number of $`s`$ quarks is conserved separately from that of the net number of light quarks. In our model, the two light quarks are both massless, and since we neglect electromagnetic interactions, there is no difference between $`u`$ and $`d`$ quarks.

### A Quark-hadron phase equilibrium

To study the phase equilibrium between the hadrons and the QGP, we implement the Gibbs conditions on the volume (or bulk) part of the thermodynamic potential, i.e. also in the QGP density of states we only use the terms proportional to the volume. (We shall see later that the surface terms give rise to interesting features near bulk phase equilibrium.) In this way, we obtain the phase diagram in Fig. 2. The QGP phase (above the phase equilibrium lines) occupies a larger and larger region of the phase diagram as the $`s`$-quark chemical potential is raised. As one might have expected, the phase equilibrium line in the case of equal chemical potentials lies somewhere between the lines of varying $`\mu _s`$.

### B The thermodynamic potential

In this section we look at the thermodynamic potential in the vicinity of equilibrium, now using the full expression of the density of states of the QGP (i.e. including the surface terms). In Figs. 36, there is a single chemical potential, $`\mu _q`$, common to all three quark flavors.

First we show Fig. 3, the thermodynamic potential of a QGP droplet in equilibrium with an extended hadron phase near the bulk equilibrium point $`(\mu _q,T)=(320,106.1)`$MeV. The picture is as expected for a first order phase transition: Below the bulk equilibrium temperature no stable QGP droplet can form, and even somewhat above this temperature there is an energy barrier for the system to pass before the true minimum at $`R=\mathrm{}`$ can be reached. The surface terms are responsible for this energy barrier. Precisely at the phase transition temperature the volume terms cancel (by definition), and the thermodynamic potential goes to infinity as $`R^2`$ (the leading surface term).

The more interesting conclusions are reached when one considers the reverse situation, namely a hadron bubble in equilibrium with an extended phase of QGP. Figs. 4, 5 and 6 show the thermodynamic potential of this configuration near different points on the phase equilibrium line. In all three cases the potential exhibits a minimum at a radius of $`R1`$ fm of approximately the same depth, even at temperatures above phase equilibrium temperature. Such a minimum would not be present in the standard textbook treatment of phase transitions, in which a first order phase transition usually is described in terms of a phenomenological free energy containing only volume- and (positive) surface terms.

Thus, a QGP in thermodynamic equilibrium apparently contains bubbles of hadrons, and the transition from QGP to hadrons will proceed (in addition to ordinary bubble nucleation and subsequent expansion of these bubbles) by smooth expansion of the pre-formed bubbles. In addition to this smooth growth as the minimum tends to larger radii when the temperature drops, pre-formed bubbles can also grow discontinuously by passing the barrier from $`R>0`$ to some $`R`$ outside the barrier. This possibility resembles very much the ordinary nucleation scenario with a modified barrier height. At the parameters of Figs. 4, 5 and 6, the mean number of quarks in a hadron bubble of radius $`R=1`$ fm is typically 3–4, suggesting that formation of a real hadron is not a very rare event. However, the minimum being of modest depth, this interesting effect could be due to the model’s inaccurate representation of QCD.

The fact that there is no energy barrier for a hadron bubble to form in QGP does not imply that the phase transition is second order, or a smooth cross-over. In the model adopted here, the transition is inherently first order, since there is an a priori difference in entropy between the two bulk phases, and thus a latent heat. As seen from e.g. Fig. 4 it is not energetically favorable for bulk hadronic matter to form at temperatures above the phase transition temperature $`T_0`$. But the formation of isolated bubbles of hadronic matter is favorable above $`T_0`$, so the transition will appear smoother than normally expected for a first order transition, with a more gradual release of latent heat.

We now investigate the effect of letting the $`s`$ quark have its own chemical potential. To this end we show Figs. 7 and 8. In Fig. 7 the $`s`$ quark chemical potential is 4 times greater than the light quark chemical potential, whereas in Fig. 8 it is converse. The conclusion to be drawn from these figures is that a large $`s`$ quark chemical potential tends to wash out the minimum in the thermodynamic potential, present when the light quark chemical potential is less than or equal to the $`s`$ quark chemical potential. This behavior can be traced back to the positive surface tension contribution from the massive $`s`$-quarks to the thermodynamic potential.

Finally, we underline that the interesting effect of bubble formation is due to the fact that the presence of a surface alters the density of states of particles. Although these finite-size effects are negligible at large radii of the bubble or droplet, they may have important implications for the phase transition as a whole.

## IV Conclusion

We have seen that the thermodynamic potential of a hadron bubble embedded in quark-gluon plasma, exhibits a minimum at a radius $`R1`$ fm, even at temperatures somewhat above the bulk transition temperature. Thus, within the model described here, a homogeneous plasma of size larger than a few fermi is an impossibility.

Regarding the relevance of these results in connection with current and forthcoming ultrarelativistic heavy ion collisions, we conclude that if indeed a quark-gluon plasma is formed in the course of such a collision, then (according to the model considered here) this plasma phase will contain bubbles of hadrons of radius $`R1`$ fm. If such hadronic bubbles form inside the plasma, the observable effects are likely to include some blurring of the plasma signatures, since this reduces the effective plasma volume.

We emphasize that these conclusions may well be model dependent, inasmuch the minimum of the thermodynamic potential is of modest depth. However, the phenomenon of spontaneous creation of stable hadronic bubbles in a quark-gluon plasma, does seem to be well established within the the model discussed here . The main question is now, whether these hadronic bubbles are of physical nature, or merely an artifact of the model. Certainly, it would improve confidence in these results if other models of QCD, and eventually lattice calculations, were to yield similar results.

###### Acknowledgements.

JM was supported in part by the Theoretical Astrophysics Center under the Danish National Research Foundation. We thank Michael Christiansen for useful discussions.