text:
# Determination of fractal dimensions of solar radio bursts

## 1 Introduction

Nonlinear time series analysis based on the theory of deterministic chaos has turned out to be a powerful tool in understanding complex dynamics from measurements and observational time series. In particular it can provide descriptions and interpretations for irregular times series, which nevertheless might not be governed by a stochastic physical process and which are only poorly understood by linear methods. A number of recent reviews and conference proceedings shows the great interest in the field of nonlinear time series analysis (see, for instance, Grassberger et al. GrassbergerEtal91 (1991); Casdagli & Eubank CasdagliEubank92 (1992); Weigend & Gershenfeld WeigendGershenfeld93 (1993); Kugiumtzis 1994a , 1994b ; Abarbanel Abarbanel96 (1996); Kantz & Schreiber KantzSchreiber97 (1997); Schreiber Schreiber99 (1999)).

Since the development of chaos theory, it is well known that even simple dynamical systems, described by few nonlinear differential equations, can reveal a complex and quasi-irregular behavior. A central concept to characterize such systems is the so-called attractor. Under the dynamics of a deterministic system the trajectories do not cover the whole phase space, but, after all transient phenomena have faded out, converge to a subset of the phase space, the attractor. The attractor itself is invariant to the dynamical evolution. Simple examples of attractors are fixed points and limit cycles. However, when the related dynamical system is chaotic, the attractor can have a complex geometry with a fractal, i.e., non-integer, dimension. Different invariant parameters exist to describe the geometry and the dynamics of an attractor, such as dimensions, Lyapunov exponents and entropies. Besides the fractal geometry, chaotic systems have the striking property that initially neighboring trajectories diverge exponentially under the dynamics, and the growth rate is given by the Lyapunov exponent. This phenomenon results from the folding and stretching of the trajectories under the dynamics, the folding leading to the convergence of the trajectories to the attractor and the stretching to the divergence in certain directions. While the average stretching rate is given by the Lyapunov exponent, the loss of information due to the folding is quantified by the entropy. Lyapunov exponents and entropies characterize the dynamics on the attractor, and the dimension characterizes its geometry. The physical meaning of the dimension of an attractor is that it corresponds to the degree of freedom of the related dynamical system, i.e., in a deterministic case, the minimum number of ordinary differential equations needed to fully describe the system. Deterministic systems are characterized by a finite dimension. Deterministic chaotic systems have the additional characteristic that their dimension is fractal. Contrary to that, a stochastic system is characterized by an infinite dimension, indicating its infinite degree of freedom. Therefore, the determination of the dimension of an attractor enables to discriminate whether a dynamical system is deterministic or stochastic.

Previous papers exist concerning the investigation of fractal dimensions of solar radio bursts. It has to be noted, that the time series used by the different authors are not directly comparable as they represent different types of radio events. Kurths & Herzel (KurthsHerzel86 (1986), KurthsHerzel87 (1987)), Kurths & Karlick$`\stackrel{´}{\mathrm{y}}`$ (KurthsKarlicky89 (1989)), and Kurths et al. (KurthsEtal91 (1991)) analyzed decimetric pulsations and ascertained finite dimension values. Contrary to that, Isliker (1992b ) and Isliker & Benz (1994a , 1994b ) investigated different types of solar radio bursts in the metric (m) and decimetric (dm) wavelength range (type I storms, type II bursts, type III bursts, type IV events, and narrowband spikes), which did not reveal any hints for low-dimensional determinism.111The reported finite dimension for one of the analyzed narrowband spike events in Isliker (1992b ) was revised in a later paper (Isliker & Benz 1994a ).

However, these investigations rely all on the correlation dimension method. The present paper additionally introduces a complementary dimension analysis, motivated by the fact that solar radio bursts represent real-world systems, which implies some major restrictions. First, the time series cannot be expected to be stationary, and second, pure determinism is rather unlikely to be realized. Therefore we do not only concentrate on the usual way of looking at the problem: “Does the analyzed time series represent a deterministic or a stochastic system?” but in particular focus the question: “What statistical description can be extracted from a dimension analysis of the time series?” With such refined formulation of the problem, we try to make use of the concepts and tools of nonlinear time series analysis even in cases in which the determination of invariants of the dynamics, as, e.g., attractor dimensions, possibly fails.

The paper is structured as follows. Sect. 2 explains the used methods and discusses critical points in the determination of fractal dimensions from time series. In Sect. 3 the investigated data sets are characterized and the analysis procedure is described. Sect. 4 presents the results of the dimension analysis, which are discussed in Sect. 5. Finally, the conclusions are drawn in Sect. 6.

## 2 Methods

### 2.1 Phase space reconstruction

Generally, not all relevant parameters of the dynamics of a system are measured during an observation but only a one-dimensional time series is given. To reconstruct the phase space of the related dynamical system, techniques have to be applied to unfold the multi-dimensional attractor from a scalar time series. By the technique of time delayed coordinates (Takens Takens81 (1981)), from a given one-dimensional time series $`\{x(t_i)\}`$ an $`m`$-dimensional phase space $`\{𝝃_i\}`$ is built up by the prescription

$$𝝃_i=\{x(t_i),x(t_i+\tau ),\mathrm{},x(t_i+(m1)\tau )\},$$
(1)
where $`\tau `$ is the time delay.

According to the embedding theorem of Takens (Takens81 (1981)), the embedding of the attractor in the $`m`$-dimensional reconstructed phase space can be ensured if $`m(2D+1)`$, with $`D`$ the dimension of the original phase space. For time series of infinite length and accuracy, $`\tau `$ can be chosen in a more or less arbitrary way without affecting the results. However, in practice not every value for the time delay $`\tau `$ will be suitable. Too small $`\tau `$ will build up coordinates which are too strongly correlated, while for large $`\tau `$ the vector components show no causal connection. The choice of the time delay $`\tau `$ in the reconstruction of the phase space usually strongly affects the quality of the analysis, and different procedures have been worked out to ensure a proper choice of $`\tau `$. The most prominent methods use the auto-correlation time ($`1/e`$ decay time, first zero crossing, first minimum) or the first minimum of the mutual information (Fraser & Swinney FraserSwinney86 (1986)).

One advantage of the mutual information over the auto-correlation function is that it takes into account nonlinear properties of the data. The mutual information is based on the Shannon entropy (Shannon & Weaver ShannonWeaver62 (1962)), and gives the information about the state of a system at time $`t+\tau `$ that we already possess if we know the state at time $`t`$. The choice of the first minimum of the mutual information for the time delay $`\tau `$ is motivated by the fact that two successive delay coordinates should be as independent as possible without making $`\tau `$ too large.

### 2.2 Correlation dimension

The correlation dimension is one out of many definitions of fractal dimensions, and was introduced by Grassberger & Procaccia (1983a , 1983b ) to determine fractal dimensions from time series. The correlation dimension is based on distance measurements of points in phase space. Therefore, as first step, from the time series $`\{x(t_i)\}`$ the phase space vectors $`\{𝝃_i\}`$ have to be constructed. With the reconstructed vectors, the correlation integral $`C(r)`$ can be calculated, which is given by the normalized number of pairs of points within a distance $`r`$. As the correlation dimension is based on spatial correlations in phase space, it is an important precaution to exclude serially correlated points in counting the pairs (for details see Sect. 2.4.2), and the length of the window, $`W`$, should at least cover all points within the auto-correlation time (Theiler Theiler86 (1986)). With this correction, the correlation integral is given by

$$C(r)=\frac{2}{\mathrm{N}_{\mathrm{pairs}}}\underset{i=1}{\overset{\mathrm{N}W}{}}\underset{j=i+W}{\overset{\mathrm{N}}{}}\mathrm{\Theta }(r𝝃_i𝝃_j),$$
(2)
$$\mathrm{N}_{\mathrm{pairs}}=(\mathrm{N}W)(\mathrm{N}W+1),$$
(3)
where N denotes the overall number of data points, and $`\mathrm{\Theta }`$ is the Heaviside step function. For small distances $`r`$, the correlation integral $`C(r)`$ is expected to scale with a power of $`r`$, and the scaling exponent defines the correlation dimension $`D_c`$:

$$C(r)r^{D_c},\mathrm{for}r0,$$
(4)
$$D_c=\underset{r0}{lim}\frac{\mathrm{ln}C(r)}{\mathrm{ln}r}.$$
(5)
Practically one computes the correlation integral for increasing embedding dimension $`m`$ and calculates the related $`D_c(m)`$ in the scaling region. If the $`D_c(m)`$ reach a saturation value $`D_c`$ for relatively small $`m`$, this gives an indication that an attractor with dimension $`D_c`$ exists underlying the analyzed time series.

### 2.3 Local pointwise dimensions

The local pointwise dimension is a locally defined variant of the correlation dimension. Its definition is based on the probability $`p_i(r)`$ to find points in a neighborhood of a point $`𝝃_i`$ with size $`r`$:

$$p_i(r)=\frac{1}{\mathrm{N}_{\mathrm{pairs}}}\underset{\genfrac{}{}{0pt}{}{j=1}{|ji|W}}{\overset{\mathrm{N}}{}}\mathrm{\Theta }(r𝝃_i𝝃_j),$$
(6)
where $`\mathrm{N}_{\mathrm{pairs}}`$ gives the actual number of pairs of points in the sum. For small distances $`r`$, $`p_i(r)`$ is expected to scale with a power of $`r`$, and the scaling exponent $`D_p(𝝃_i)`$ gives the local pointwise dimension at point $`𝝃_i`$:

$$p_i(r)r^{D_p(𝝃_i)},\mathrm{for}r0,$$
(7)
$$D_p(𝝃_i)=\underset{r0}{lim}\frac{\mathrm{ln}p_i(r)}{\mathrm{ln}r}.$$
(8)
Averaging $`D_p(𝝃_i)`$ over all points of the time series or a number of reference points yields the averaged pointwise dimension, $`\overline{D}_p`$, which is equivalent to the correlation dimension and gives a global description of the geometry of an attractor:

$$\overline{D}_p=\frac{1}{\mathrm{N}_{\mathrm{ok}}}\underset{i=1}{\overset{\mathrm{N}_{\mathrm{ok}}}{}}D_p(𝝃_i),$$
(9)
with $`\mathrm{N}_{\mathrm{ok}}`$ the number of accepted reference points.

Since the $`D_p(𝝃_i)`$ are local functions and defined for each point, they are, on the one hand, a function of the position $`𝝃_i`$ on the attractor, characterizing its local geometry. On the other hand, the local dimensions can be interpreted as a function of time, $`D_p(t_i)`$, since they reflect the temporal evolution, in which the points $`𝝃_i`$ on the attractor are covered by the dynamics (Mayer-Kress Mayer-Kress94 (1994)). Based on this fact, local dimension estimations have the interesting property that they enable to cope with non-stationary data.

An example is shown in Fig. 1, which illustrates the local pointwise dimensions calculated for a simulated time series, made up of three sections related to different attractors, the Lorenz, a limit cycle and the Rössler attractor. The evolution of the local pointwise dimensions $`D_p(t_i)`$ detects the different attractors successively operating in time. The averaged pointwise dimensions $`\overline{D}_p`$ in the respective sections deviate less than 10% from the true values. As for the calculation of the local dimension at point $`𝝃_i`$ the distances to all points of the time series, even those related to a different attractor, are taken into account, it is a quite striking feature that the different attractors can be disentangled by the method. Moreover, for this exemplary analysis the length of the time series was taken rather short in order to mimic the conditions of observed time series. Such a behavior reveals that the reference point $`𝝃_i`$, which is equivalent to a time $`t_i`$ of the dynamical evolution of the system, dominates the local dimension calculation. However, we want to stress that small changes of the attractor dimension, as, e.g., in the case of the Lorenz and the Rössler attractor (see Fig. 1), cannot be detected .

### 2.4 Pitfalls in dimension estimations

The determination of fractal dimensions from time series, characterized by finite length and accuracy, includes many pitfalls, which can lead to quite spurious results. In this chapter the most prominent problems will be discussed and the strategies used by the authors to avoid such pitfalls. A brief but quite dense review with respect to critical points in the determination of the correlation dimension can be found in Grassberger et al. (GrassbergerEtal91 (1991)). As the local pointwise dimensions are a variant of the correlation dimension, most of the problems occur in a similar way.

#### 2.4.1 Noise and quality of the scaling region

As expressed in Eq. 4 for the correlation integral $`C(r)`$ and Eq. 7 for the probability $`p_i(r)`$, respectively, a scaling behavior is expected for $`r0`$. However, for real time series, which are contaminated by noise and which are of finite length, the scaling behavior is expected to occur at intermediate length scales. Noise is acting at small scales and therefore dominating the scaling behavior for small $`r`$. Deviations from ideal scaling at large length scales are due to edge effects caused by the finite length of the time series.

Fig. 2 shows such typical scaling behavior for a time series of finite length and accuracy. For better illustration we have plotted the local slopes of the correlation integral, given by the expression

$$\nu (r)=\frac{\mathrm{d}\mathrm{ln}C(r)}{\mathrm{d}\mathrm{ln}r}.$$
(10)
In the case of ideal scaling, the $`\nu (r)`$, calculated for different embedding dimensions $`m`$, form straight lines parallel to the $`x`$-axes, the so-called plateau region, with the constant $`y`$-value corresponding to $`D_c(m)`$. However, it is typical for observational time series that at least three different parts in the $`\nu (r)`$ curves are distinguishable, described in the caption of Fig. 2. To avoid spurious finite dimensions which might arise from a misinterpretation of the existence or non-existence of a physically relevant scaling region, we implemented an algorithm to automatically check for scaling behavior (see Sect. 2.5).

#### 2.4.2 Temporal correlations

The correlation integral $`C(r)`$ and the probability $`p_i(r)`$ are a measure of spatial correlations on the attractor. They are basically calculated by counting pairs of points which are closer to each other than a given distance $`r`$. However, at small scales successive points of the time series give an additional contribution, since they are close in time. Nevertheless, such points do not reflect spatial correlations, i.e., a clustering of points in phase space, but are serially correlated by the temporal order of the time series. To avoid the spurious contribution of serially correlated points, which can cause strong artificial effects, especially for data recorded with a high sampling rate, at least all pairs of points closer than the auto-correlation time have to be excluded (Theiler Theiler86 (1986)):

$$W>t_{corr}.$$
(11)
A similar phenomenon can occur when the analyzed time series is rather short. If the time series is too short to ensure that the attractor is well covered with points, then most of the points of the series are serially correlated, which again might result in spurious dimensions. Different relations have been derived for the minimum length of a time series needed for dimension estimations, as the one by Eckmann & Ruelle (EckmannRuelle92 (1992)):

$$D_c<2\mathrm{log}_{10}\mathrm{N},$$
(12)
with N the length of the time series. However, not only the length of the data series is of relevance but also the length with regard to the sampling rate. Applications to different kind of time series, from known chaotic attractors as well as measured time series (Brandstater & Swinney BrandstaterSwinney87 (1987); Kurths et al. KurthsEtal91 (1991); Isliker 1992a ; Isliker & Benz 1994a ), have revealed that the analyzed time series should at least cover 50 structures – “structure” meaning a full orbit in phase space or generally a typical time scale of the analyzed time series. According to Isliker (1992a ), we define the structures by the first minimum of the auto-correlation function,

$$\mathrm{N}_{str}50\mathrm{with}\mathrm{N}_{str}:=\frac{\mathrm{N}\mathrm{\Delta }}{t_{corr}},$$
(13)
where $`\mathrm{N}_{str}`$ gives the number of structures, N the length of the time series, $`\mathrm{\Delta }`$ the temporal resolution, and $`t_{corr}`$ the auto-correlation time.

#### 2.4.3 Stationarity

As shown by Osborne et al. (OsborneEtal86 (1986)) and Osborne & Provenzale (OsborneProvenzale89 (1989)), the determination of the correlation dimension from non-stationary stochastic processes, can erroneously lead to finite dimension values. To take into account that problem, we applied a stationarity test proposed by Isliker & Kurths (IslikerKurths93 (1993)), which is based on the binned probability distribution of the data. To check for stationarity one divides the time series into subsections, and compares the probability distribution of the section under investigation with the probability distribution of the first half of it by a $`\chi ^2`$-test. By the use of this test we searched for stationary subsections in the radio burst time series, and only to such stationary subsections the correlation dimension analysis was applied.

However, with the concept of local dimensions it is possible to cope with non-stationary data. One important advantage of the local dimension method is that it can enable to detect dynamical changes in a time series. To make use of this potentiality, we calculated the local pointwise dimensions from the whole time series instead of using stationary subsections. Moreover, since the statistics in the calculation of the local pointwise dimensions grows linearly with the length of the time series, whereas the correlation dimension grows with the square of the number of points, for the local pointwise dimension analysis the time series should be kept as long as possible. However, to avoid spurious results due to non-stationarities, we applied a surrogate data test (Sect. 2.6).

#### 2.4.4 Intermittency

Intermittency describes the phenomenon that a time series is interrupted by quiet phases or phases of very low amplitudes. Such a phenomenon is quite typical for chaotic systems, but can cause problematic situations when calculating fractal dimensions from limited time series. In phase space the intermittent sections represent regions, which, in the context of the global attractor scale, degenerate to a point. This would trivially result in an erroneously low dimension value. One way to cope with this problem is just to discard intermittent phases from the analyzed time series. Since the dimension of an attractor is a geometric descriptor, which means a static quantity, it is not influenced by the serial order of points, and therefore discarding subsections is a valid strategy.

As an example, the top panel of Fig. 3 shows the subsection of a type I burst series. The middle and bottom panels show the results of the correlation dimension analysis, calculated from the whole time series (middle panels) and the time series after the intermittent section was eliminated (bottom panels). A comparison of the middle and bottom panels reveals that the intermittent phase causes a deformation in the curves of the correlation integral and the related local slopes. Eliminating the intermittent section, these deformations disappear. However, we want to stress that the used algorithm for the automatic detection of the scaling region does not identify the sink-region as a scaling region, which could result in spurious finite dimensions. Therefore, by the set up of the scaling algorithm we avoid erroneous low dimensions caused by intermittency effects.

### 2.5 Scaling and convergence test

The practical computation of the local pointwise dimensions turns out to be more difficult than the correlation dimension, particularly as for each reference point the scaling region has to be determined and the related $`D_p(𝝃_i,m)`$, which are calculated in the scaling region, have to be checked for convergence with increasing embedding dimension $`m`$. For this purpose, an automatic and fast procedure is needed. Moreover, such an automatic procedure can also be used in the correlation dimension analysis to avoid subjective influences on the scaling and convergence judgement. We implemented such an algorithm, based on the one used by Skinner et al. (SkinnerEtal91 (1991)), but modified in order to reach higher stability and significance. This automatic procedure searches for the scaling region, defined as the longest linear range in the $`\mathrm{ln}p_i(r)`$ curves, tests if the scaling region is of significant length, and finally checks if the $`D_p(𝝃_i,m)`$ are converging with increasing $`m`$. Only for points $`𝝃_i`$, which pass the scaling and the convergence test, a local pointwise dimension $`D_p(𝝃_i)`$ is accepted. In the following we give a description of the algorithm.

According to Eq. 6, for each reference point $`𝝃_i`$ and each considered embedding dimension $`m`$, we calculate the cumulative histogram of the $`p_i(r)`$, dividing the $`r`$-range into 100 equidistant points in the logarithmic representation, which cover the whole range from the smallest to the largest actual distance. The local slopes of the $`\mathrm{ln}p_i(r)`$ versus $`\mathrm{ln}(r)`$ curves, given by

$$\nu _i(r)=\frac{\mathrm{d}\mathrm{ln}p_i(r)}{\mathrm{d}\mathrm{ln}r},$$
(14)
calculated in the scaling region, represent the $`D_p(𝝃_i,m)`$. For the determination of the location of the scaling range, we shift a window with length 10 points through the overall $`r`$-range, and for each of these windows a least squares linear fit is applied to the $`\mathrm{ln}p_i(r)`$ versus $`\mathrm{ln}r`$ curves. The slope of the linear fit corresponds to the average value of the local slopes in the corresponding $`r`$-range, denoted as $`\overline{\nu }_i(r)`$. Starting with the first window, successive $`\overline{\nu }_i(r)`$ are compared with the corresponding quantity in the first window, until the difference is larger than a certain threshold value, chosen as 20% of the initial $`\overline{\nu }_i(r)`$. In this case, the position of the first point of the start window and the last point of the actual window are stored, and the procedure continues with the second window as start window. If the new positions determined cover a larger range than the old ones, the old values are overwritten by the new ones, and so on. The two stored values remaining at the end give the location of the scaling region, in which we calculate $`D_p(𝝃_i,m)`$. If the determined scaling region is smaller than 20% of the overall length, it is interpreted as not significant and rejected. To avoid spurious results at large $`r`$, which usually correspond to small $`\overline{\nu }_i(r)`$ values, in the whole procedure we suppress values $`\overline{\nu }_i(r)<2`$. Fig. 4 shows a sample application of the algorithm detecting the scaling range.

After the scaling region is determined for each point $`𝝃_i`$ and each embedding dimension $`m`$, for those points which reveal a scaling region of significant length for all considered $`m`$-values, we test the convergence of the $`D_p(𝝃_i,m)`$ with increasing $`m`$. This is simply done by averaging the $`D_p(𝝃_i,m)`$ over four successive $`m`$-values. If the standard deviation turns out to be less than 15% of the average value, the convergence is accepted, and the average of the $`D_p(𝝃_i,m)`$ over the considered $`m`$-range is taken as local pointwise dimension $`D_p(𝝃_i)`$. The different kind of threshold values used have been adjusted by application of the algorithm to different time series, and slight changes do not qualitatively change the outcome.

### 2.6 Surrogate data test

As already mentioned, local dimensions enable to deal with non-stationary data. However, to avoid spurious dimensions which can result from non-stationary stochastic processes, we applied a surrogate data test (Osborne et al. OsborneEtal86 (1986); Theiler et al. TheilerEtal92 (1992)). For this purpose, a Fourier transform of the time series is performed, the Fourier phases are randomized, and finally to this phase-altered series the inverse Fourier tranform is applied. The phase randomization keeps the power spectrum unchanged while linear correlations in the time series are eliminated. If the results of the dimension analysis of such an ensemble of surrogate data are significantly different from those computed from the original data, the hypothesis can be rejected that the obtained results are caused by a linear stochastic process, and nonlinearity in the data is detected.

To quantify the statistical significance we make use of a hypothesis testing given in Theiler et al. (TheilerEtal92 (1992)). The null hypothesis assumes that the original data represent linearly correlated noise. To compute a discrimination statistics we need a set of surrogate data, from which the same statistical quantities are derived as from the original time series. In particular, we use the averaged pointwise dimensions as function of the embedding dimension, $`\overline{D}_p(m)`$. Let $`Q_o`$ denote the statistical quantity computed for the original time series, and $`Q_{s_i}`$ for the $`i`$th surrogate generated under the null hypothesis. Let $`\mu _s`$ and $`\sigma _s`$ represent the (sample) mean and standard deviation of the distribution of $`Q_s`$. With this notation the measure of “significance”, $`S`$, is given by the expression

$$S=\frac{|Q_o\mu _s|}{\sigma _s}.$$
(15)
## 3 Data sets and analysis procedure

The data sets are single frequency recordings from the multichannel radio-polarimeter of the Trieste Astronomical Observatory, which is operating in the dm-m wavelength range. The investigated data sets are recorded at the frequencies 237, 327, 408 and 610 MHz, with a sampling rate of 50 Hz, i.e. a temporal resolution $`\mathrm{\Delta }=20`$ ms. We analyzed 30 data sets of type I storms and 27 sets of type IV events, which cover samples with different kind of fine structures, such as pulsations, fast pulsations, sudden reductions and spikes. Some of the events were analyzed at different times and/or frequencies. Therefore, the 30 data sets of type I storms represent 24 different events, the 27 data sets of type IV events 20 different ones. In Table 1 and Table 2 we give a description of the type I and type IV events, respectively. The main criteria for the data selection from the solar radio burst data archive of the Trieste Observatory were:

1. The selected data sets were representative for the particular types of events.
2. To ensure a high signal-to-noise ratio only intense events were selected.
3. The related time series were substantially long and fulfilled Eq. 13.

For the analysis, the predominant polarization sense, LCP (Left-handed Circular Polarization) or RCP (Right-handed Circular Polarization) of the burst series was used.

The first step in the correlation dimension analysis was to search for stationary subsections by shifting windows with decreasing length through the time series and applying the stationarity test proposed by Isliker & Kurths (IslikerKurths93 (1993)). Only those stationary subsections which still fulfilled the minimum length criterion of Eq. 13 were accepted for further analysis, and the correlation dimension was calculated only from such subsections. The analysis was repeated with different values for the delay parameter$`\tau `$, located around the first minimum of the mutual information. The relevant quantities were calculated up to embedding dimension $`m=20`$. Finally, the algorithm for automatically searching the scaling region, checking its validity and testing the convergence behavior was applied. The convergence was checked for $`m`$-intervals containing four successive embedding dimensions, starting with $`m=10`$.

The pointwise dimension analysis was basically carried out in the same manner, except that the overall time series was used instead of stationary subsections. Moreover, only for points which passed both scaling and convergence test a local pointwise dimension was accepted. Finally, the pointwise dimension analysis was repeated for 10 different sets of surrogate data to test against the null hypothesis that the results are caused by linearly correlated noise, and to get evidence on nonlinearity in the data.

## 4 Results

### 4.1 Correlation dimension analysis

The application of the stationarity test led to the result, that only two of the preselected events did not reveal stationary subsections fulfilling Eq. 13. These events were excluded from further analysis. For each of the other events, the longest stationary subsections analyzed are listed in Table 1 and Table 2. If more stationary subsections were found, for the analysis the longest three were selected, and non overlapping or only partially overlapping subsections were preferred. The correlation dimension analysis of the stationary subsections did not reveal a finite dimension for any of the events. We identify mainly three cases, in which negative results occurred: no convergence, no scaling region, and deformed scaling region. These three cases are of course ideal ones, and the practice quite often revealed a mixture of it. Therefore, for the single event sections analyzed, we do not specify the different reasons causing the negative result, but describe and discuss the principal cases in the following subsections in a general frame.

#### 4.1.1 No convergence

If the $`D_c(m)`$ converge with increasing embedding dimension $`m`$, this gives indications for the existence of a low-dimensional attractor underlying a time series. In the analyzed data sets, for which the $`C(r)`$ curves revealed a clear scaling region, no convergence with increasing $`m`$ occurred. As an example, Fig. 5 shows the divergent behavior for a type IV event. Even for rather high embedding dimensions a distinct scaling region exists. However, as the curves of the local slopes $`\nu (r)`$ clearly reveal, the plateau of the scaling region moves to higher values for increasing $`m`$, indicating that the $`D_c(m)`$ are divergent. The meaning of such a divergent behavior can be manifold, being related to the physical state of the system as well as to restrictions regarding the analysis methodology:

1. The underlying physical system is stochastic.
2. The signal is the output of a deterministic but high-dimensional system, with a dimension too high to be extracted from the given time series of finite length.
3. The signal represents a system to which noise is coupled intrinsically to the dynamics.
4. The analyzed time series is the result of different physical systems which are independently operating at the same time, e.g., different uncoupled radio burst sources simultaneously present on the sun, whose emissions sum up to the measured signal.
5. Despite the careful data selection, the measurement noise is still too high for the kind of analysis carried out and dominates the results.
6. The choice of the time delay $`\tau `$, which is a quite critical and sensitive parameter in the correlation dimension analysis, is not optimal.

#### 4.1.2 No scaling region

There exist cases, in which for increasing embedding dimension $`m`$ the scaling region disappears. In the curves of the correlation integral $`C(r)`$ no linear range can be detected and, correspondingly, in the graphs of the local slopes $`\nu (r)`$ no plateau region exists. In such cases the quantity $`D_c(m)`$ is not defined. An example is shown in Fig. 6. The main reason for such a behavior might be that the analyzed data set is too short and/or too strongly contaminated by noise. In such a case the deviations from ideal scaling at small and large length scales $`r`$ merge at intermediate $`r`$, and cause the scaling region to disappear. This effect is stronger for higher embedding dimensions, since for increasing $`m`$ the attractor is sparsely covered with points, which can be verified in Fig. 6. As shown by Schreiber & Kantz (SchreiberKantz95 (1995)), even small amounts of measurement noise can conceal possible scaling behavior.

#### 4.1.3 Deformed scaling region

The third case, in which the correlation dimension analysis led to negative results, is, that the curves of the correlation integral are highly deformed. Such deviations from ideal scaling are mainly caused by the presence of intermittent sections, or generally, by the presence of very different amplitude scales in the data. Fig. 3 shows such deformed scaling regions for a stationary subsection of a sample type I storm, caused by an intermittent section. This can be clarified as the elimination of the intermittent section causes the deformation to disappear. Explained by the typical presence of intermittent phases in type I storms, the correlation dimension analysis rather often results in deformed plateau regions for type I burst series.

### 4.2 Local pointwise dimension analysis

The third set of columns in Table 1 and Table 2 contains the results of the pointwise dimension analysis for the type I and the type IV events, respectively. $`\overline{D}_p`$ gives the averaged pointwise dimension with standard deviation $`\sigma `$, calculated over an embedding range $`m=1013`$, and “ok” denotes the percentage of points which passed the scaling and convergence test. $`\mathrm{\Delta }\overline{D}_p`$ gives the increase of the averaged pointwise dimension, when the embedding range is increased by two dimensions, from $`m=1013`$ to $`m=1215`$. “surr.” denotes the outcome of the surrogate data test, and positive means that the null hypothesis, which states that the results are caused by linearly correlated noise, can be rejected, and that significant evidence for nonlinearity in the data is given.

As an example, Fig. 7 illustrates the results of the local pointwise dimension analysis obtained for a type IV event with quasi-periodic pulsations. The top panel shows the time series, the middle panel the time evolution of the local pointwise dimensions, $`D_p(t_i)`$. In the bottom panels the histograms of the local dimensions are plotted, calculated over different embedding ranges. About half of the points of the time series passed the scaling and convergence test and were used to compute the averaged pointwise dimension. However, as the histograms of the local dimensions calculated over increasing embedding dimensions ($`m=710`$, $`m=912`$, $`m=1114`$) reveal, no absolute convergence exists, but a slight increase with increasing $`m`$-ranges occurs. This can be clearly seen in the right bottom panel of Fig. 7, as for higher embedding ranges the center of the histogram moves to higher dimension values. The same phenomenon occurs for all analyzed samples, as the positive $`\mathrm{\Delta }\overline{D}_p`$ values indicate, on the average $`\mathrm{\Delta }\overline{D}_p5\%`$.

This phenomenon, on the one hand, could result from the fact that the $`D_p(𝝃_i)`$ do not really converge with increasing $`m`$ and no convergence to a low-dimensional attractor exists. On the other hand, a comparison with simulated time series from well known chaotic attractors contaminated with Gaussian noise reveals a similar behavior. For selected samples we repeated the analysis after application of a simple nonlinear noise reduction to the data (Schreiber Schreiber93 (1993)), to find out if the increase is caused by measurement noise. As it is common in nonlinear time series analysis, the used noise reduction method does not rely on frequency information but makes use of the structure in the reconstructed phase space. Fig. 8 shows an example of the application of the noise reduction method. Repeating the analysis with the times series after the application of the noise reduction led to the effect, that the dimension increase is softened but not fully eliminated though. Such result suggests that the dimension increase is not caused by measurement noise contaminating the signal of a deterministic system.

The surrogate data analysis yields a positive outcome for most of the analyzed samples, indicating that the obtained results are not an artifact, which could arise when the analysis is applied to non-stationary stochastic data. For $``$ 75% of the type I and $``$ 85% of the type IV burst time series with a significance of $`3\sigma `$, we can reject the null hypothesis that the results are caused by a linear stochastic process, and have evidence for nonlinearity in the data. Fig. 9 illustrates the outcome of the surrogate data analysis for a sample type I storm. The top panel shows the averaged pointwise dimension as a function of the embedding dimension $`m`$, calculated from all points passing the scaling test. Diamonds represent the original data, crosses the surrogates. The bottom panel shows the related significance $`S`$, which reaches at the maximum a level of $`20\sigma `$, giving strong evidence for nonlinearity in the data. However, corresponding to the positive $`\mathrm{\Delta }\overline{D}_p`$ values, the $`\overline{D}_p(m)`$ of the original data slightly increase with increasing $`m`$ and do not show a definitive convergence to a finite dimension value.

## 5 Discussion

The outcome of the dimension analysis does not allow to claim low-dimensional determinism for the analyzed data sets. First, the correlation dimension analysis failed in all cases. Second, the obtained averaged pointwise dimension values $`\overline{D}_p`$ are too high to characterize a low-dimensional physical system. Third, the $`\overline{D}_p`$ do not reveal an absolute convergence with increasing embedding dimension $`m`$. On the other hand though, the local pointwise dimensions obtained over a specific embedding range yield a quite distinct behavior, and the surrogate data analysis evidences nonlinearity in the data.

In Fig. 10 we show a comparison of the correlation dimension and the local pointwise dimension analysis of a sample type IV event. The correlation integral and the related local slopes do not give any evidence for low-dimensional determinism, since no significant scaling region and no convergent behavior occurs. However, the corresponding curves of the local pointwise dimension analysis reveal a distinct scaling region and a clear convergence to a finite and low dimension value for certain points $`𝝃_i`$. The figure clearly illustrates that the calculated local dimensions do not represent an artifact, which might arise, for instance, if the automatic scaling and convergence procedure is not well adapted to the analysis. Moreover, such comparison suggests that the local dimension analysis of a times series is more robust than the classical correlation dimension method. The main reason might be that in the correlation dimension analysis the scaling behavior itself is a global property, since all pairs of points contribute to the correlation integral, and the scaling region can possibly be smeared out by such an averaging process. Contrary to that, the local pointwise dimensions are based on the scaling behavior at each single point, and can reveal a well defined scaling at certain points.

In Table 3 we give a summary of the averaged pointwise dimensions for the different event types. For the type IV events we additionally calculated the quantities separately for the different subtypes. We list the event type, the number of events belonging to each type (or subtype), the number of events passing the surrogate test, and the average of the pointwise dimensions over the respective event types, calculated only from the samples for which the surrogate data test gave a positive result. On the average, the type I storms reveal lower $`\overline{D}_p`$ values than the type IV events and a significant smaller standard deviation, indicating that type I storms represent a comparatively homogeneous class. The large standard deviation for the type IV burst series basically reflects the different subtypes. The average of the particular type IV subtypes yields $`\overline{D}_p`$ values which are significantly different: pulsations and sudden reductions reveal the lowest values, followed by the type IVs with fine structures of no particular kind and spikes, whereas the fast pulsations are characterized by the highest values.

Such statistics suggests that the $`\overline{D}_p`$ are quite well representing the different types of radio bursts under investigation. This means that even if low-dimensional determinism cannot be proved, the local dimension analysis can provide a quantitative description of the events, which is not possible with the commonly used correlation dimension. As argued by Schreiber (Schreiber99 (1999)), such description has the drawback that it does not provide an invariant characterization of a system. On the other hand, it offers an alternative statistical approach for systems, for which pure determinism cannot be established. This is of particular interest in astrophysics, since astrophysical systems represent real-world systems, which cannot be influenced by the observer and are highly interconnected with their surroundings, making pure determinism rather improbable.

Moreover, we claim that on a comparative basis the retrieved dimension values are related to the degree of freedom of the system. From numerical experiments with known chaotic attractors contaminated with Gaussian noise we retrieved pointwise dimensions, which are slowly increasing with increasing embedding dimension $`m`$, i.e. the absolute convergence to the definite attractor dimension disappeared due to the contamination of the deterministic signal with a stochastic component. This behavior is similar to the one of the analyzed radio events. However, from this similarity we cannot conclude that the pointwise dimension analysis is indicative for hidden deterministic chaos in the radio burst time series, since the retrieved dimension values are too high to characterize low-dimensional determinism, and for a reliable dimension analysis of high-dimensional systems much longer time series are needed than we have at disposal (– a limitation which is intrinsic to time series representing real-world systems). Nevertheless, what we can infer from this similarity is that the retrieved dimension values, even if not representing attractor dimensions, are still indicative for the degree of freedom of the physical system underlying the time series, characterizing its complexity. Based on this fact, we make use of the retrieved dimension values to describe the complexity of the related systems in a comparative way, without claiming or supposing the presence of low-dimensional determinism.

## 6 Conclusions

In the following items we give a summary of the main results obtained by the presented dimension analysis of several types of solar radio events, based on the correlation dimension and the local pointwise dimension method. The results are relevant concerning the physics of the analyzed events as well as the different methods applied.

1. The analysis does not enable to claim low-dimensional determinism in the time series. This outcome is in agreement with the results obtained by Isliker (1992b ) and Isliker & Benz (1994a , 1994b ), who also, among others, investigated type I storms, type IV events, and spikes. We cannot confirm the results of Kurths & Herzel (KurthsHerzel86 (1986), KurthsHerzel87 (1987)), Kurths & Karlick$`\stackrel{´}{\mathrm{y}}`$ (KurthsKarlicky89 (1989)), and Kurths et al. (KurthsEtal91 (1991)), who obtained finite dimension values for decimetric pulsations. However, the outcome of the present paper does not exclude deterministic chaos in the analyzed time series but makes pure low-dimensional determinism, characterized by few free parameters, rather improbable.
2. The analyzed time series are not fully stochastic, i.e. white noise. This fact we infer from the distinctly slower increase of the dimension values with increasing embedding dimension than expected for fully stochastic processes, which always fill the whole phase space, i.e. $`D(m)m`$. Moreover, the surrogate data analysis suggests that the time series do not represent linear stochastic processes.
3. For most of the analyzed data sets we have evidence that nonlinearity in the time series is present (given on a $`3\sigma `$ level by means of a surrogate data test).
4. A comparison of the two different methods used for the determination of fractal dimensions reveals that the local dimension method is more stable and enables more physical insight than the classical correlation dimension method. The local dimension analysis can provide a statistically significant quantity for systems, which cannot be characterized by invariants of the dynamics, probably since they are in fact not purely deterministic. Such quantities can be of special interest for comparative studies, investigating interrelations between different time series (which, e.g., can be useful for classificational purposes) or investigating intrarelations in between one time series (in order to detect dynamical changes).
5. The retrieved pointwise dimension values can be interpreted in terms of complexity of the underlying physical system. In this frame our analysis indicates that spikes and fast pulsations are the signature of systems of higher complexity than pulsations, sudden reductions and type I storms.

In relation with other kind of analysis of solar radio bursts the presented results might give further ideas on the physics of the events. In the following we present a short discussion in this respect, applied to pulsation and spike events, which are quite striking features associated with solar flares.

Spikes have been intensively studied during recent times. Their short duration and small bandwidth gives rise to the evidence that they are associated with the energy fragmentation process in solar flares (Benz Benz85 (1985), Benz86 (1986)). Based on this connection, Schwarz et al. (SchwarzEtAl93 (1993)) performed a nonlinear analysis by means of symbolic dynamics methods, interpreting the spikes appearance in the frequency-time domain as spatio-temporal patterns. This analysis gives indications that the simultaneous appearance of spikes at different frequencies is not a purely stochastic phenomenon but may be caused by a nonlinear deterministic (not necessarily low-dimensional) system or by a Markov process, compatible with a scenario in which spikes at nearby locations are simultaneously triggered by a common exciter, i.e. the localized sources are causally connected. In the present paper we find evidence for the spike events analyzed, that they do not represent a purely stochastic phenomenon in their temporal order either, even if the degree of freedom of the related physical system is expected to be quite high. Interpreting this result in the frame of the scenario suggested by Schwarz et al. (1993), it might give indications that the triggering of successive spikes by a localized source is not caused by a fully stochastic process, but reveals some (possibly weak) kind of nonlinear causal connection. However, this inference is restricted to the assumption that the spikes time series rather reflect the physical conditions of the triggering mechanism than those of the emission.

Pulsations, although a rather marginal phenomenon in the course of solar flares, have reached a wealth of attention, especially due to the very regular features they sometimes reveal (for a review see Aschwanden Aschwanden87 (1987)). In previous investigations of the dimensionality of solar pulsations (Kurths & Herzel KurthsHerzel86 (1986), KurthsHerzel87 (1987); Kurths & Karlick$`\stackrel{´}{\mathrm{y}}`$ KurthsKarlicky89 (1989); Kurths et al. KurthsEtal91 (1991)) the presence of low-dimensional determinism is reported, with dimensions $`2.5D3.5`$. Moreover, for one single event a dynamical evolution from a limit cycle to a low-dimensional chaotic behavior was found (Kurths & Karlick$`\stackrel{´}{\mathrm{y}}`$ KurthsKarlicky89 (1989)). Although we cannot confirm these results, we want to stress that our analysis suggests that pulsation events, especially quasi-periodic pulsations, represent the least complex phenomena among the analyzed types of radio events, i.e. their degree of freedom is expected to be lower than that of other burst types, even if not low-dimensional. The inferred high-dimensionality and the nonlinear structures detected do not match with linear MHD oscillation models for pulsations (e.g., Rosenberg Rosenberg70 (1970); Roberts et al. RobertsEtAl84 (1984)), in which only a few eigenmodes are excited, but rather favor models of self-organizing systems of plasma instabilities, which comprise periodic as well as low- and high-dimensional chaotic behavior. Such a self-organizing model for the electron-cyclotron maser instability, based on a Lotka-Volterra type equation system, is discussed in Aschwanden & Benz (AschwandenBenz88 (1988)), however restricted to limit cycle solutions.

###### Acknowledgements.

We thank P. Zlobec for helpful discussions and his support in the data classification. A.V. and A.H. acknowledge the Austrian Ministry of Sciences. M.M. acknowledges the financial support by ASI and MURST.


format: text

text:
# Depinning of elastic manifolds

## Abstract

We compute roughness exponents of elastic $`d`$–dimensional manifolds in ($`d+1`$)–dimensional embedding spaces at the depinning transition for $`d=1,\mathrm{},4`$. Our numerical method is rigorously based on a Hamiltonian formulation; it allows to determine the critical manifold in finite samples for an arbitrary convex elastic energy. For a harmonic elastic energy ($`\mathrm{\Delta }^2`$–model), we find values of the roughness exponent between the one-loop and the two-loop functional renormalization group result, in good agreement with earlier cellular automata simulations. We find that the $`\mathrm{\Delta }^2`$–model is unstable with respect both to slight stiffening and to weakening of the elastic potential. Anharmonic corrections to the elastic energy allow us to obtain the critical exponents of the quenched KPZ class.

Elastic manifolds in random media are an important issue of current research in statistical physics.Kardar ; Barabasi  In the zero-temperature motion of these manifolds, subject to a driving force $`f`$, the “depinning threshold” $`f_c`$ plays a central role: For forces $`f>f_c`$, the elastic manifold moves with finite velocity, while it is pinned for $`ff_c`$. Among the subjects studied at and around $`f_c`$ are the “creep” motion for $`f<f_c`$ at finite temperature, the scaling of the velocity for small forces $`f\stackrel{>}{}f_c`$, and the statistical properties of the pinned critical manifold at $`f_c`$, especially its roughness exponent $`\zeta `$.

If one neglects velocity-dependent terms in the equations of motion of the manifold, which one assumes to be a single-valued function $`h(𝐱,t)`$ (no overhangs), its dynamics is governed by a functional $`E(\{h,𝐱\})`$ incorporating potential energy due to the driving force $`f`$, the disorder $`\eta (𝐱,h)`$, as well as its internal elastic energy $`E_{\text{el}}`$,

$$_th(𝐱,t)=\frac{E}{h(𝐱)}=f+\eta (𝐱,h)\frac{E_{\text{el}}}{h(𝐱)}.$$
(1)
Note that $`𝐱`$ is a $`d`$–dimensional vector, in an embedding space of dimension $`d+1`$. Dimensional analysis suggests that the harmonic approximation for $`E_{\text{el}}`$ provides the only relevant term for the interplay between disorder and elasticity. This yields the quenched Edwards-Wilkinson (EW) equation

$$_th(𝐱,t)=f+\eta (𝐱,h)+a^2h(𝐱).$$
(2)
The model described by eq. (2) has been an important testing ground for the concepts and techniques originally developed in the field of critical phenomena. Functional renormalization group techniquesNattermann2 ; Narayan  were used to perturbatively compute critical exponents of the quenched EW equation, which were believed to describe generic driven manifolds at the depinning threshold.

According to the framework provided by the functional renormalization group, the manifold is flat ($`\zeta =0`$) for $`dd_{\text{uc}}=4`$. Below this upper critical dimension, the roughness exponent is expressed in an $`ϵ`$-expansion in $`ϵ=4d`$. The first order (one-loop) term Nattermann2 ; Narayan  of the expansion gives $`\zeta =ϵ/3`$. Initially, this result was believed exactNarayan  for all $`d=1,2,3`$. However, Chauve et al. Chauve2  have obtained the two-loop corrections, and found them to be non-zero.

In this work, we compute the roughness exponent of critical manifolds in finite ($`d+1`$)–dimensional samples for $`d=1,\mathrm{},4`$. Our approach is rigorously based on a Hamiltonian formulation, as we study directly the dynamics of eq. (1) for a general convex elastic energy $`E_{\text{el}}`$: A powerful numerical algorithm allows to solve this problem for a slightly modified version of eq. (1), with a discretized vector $`𝐱`$ ($`𝐱𝐢`$) and a continuous variable $`h`$. The short-range elastic energy depends on the next-neighbor distances $`\mathrm{\Delta }_{𝐢,𝜹}`$ between lattice point $`𝐢`$ and its $`2d`$ neighbors $`𝐢\pm 𝜹_j`$, as shown in Fig. 1:

$$\frac{1}{2}\underset{j=1}{\overset{d}{}}\underset{\pm 𝜹_j}{}E_{\text{el}}(|h_𝐢h_{𝐢+𝜹_j}|).$$
(3)
We introduce periodic boundary conditions for $`𝐱`$ and $`h`$ (cf Fig. 1).f-n-winding  Continuous, periodic random potentials are constructed from $`M`$ normally distributed random variables.RossoKrauth3 ; RossoKrauth5  In such a system, the sample-dependent critical force $`f_c`$ is well defined: pinned configuration with $`v_𝐢=0𝐢`$ do (do not) exist for driving forces smaller (larger) than $`f_c`$. Our algorithmRossoKrauth3  allows to decide quickly whether pinned configurations exist at a given driving force, and then zooms in on the critical force and the critical manifold.

For the first time, we are thus able to unambiguously compute the object of most direct theoretical interest in this problem, namely the critical manifold, in dimensions $`d=2,3,4`$. We computeRossoKrauth2  the roughness exponent $`\zeta `$ from the disorder-averaged mean square deviations of the critical manifolds $`h^c`$:

$$W^2(L)=\overline{(h^ch^c)^2}L^{2\zeta }\text{for}L\mathrm{}.$$
(4)
As the width of the manifold is of order $`L^\zeta `$, we scale the lateral extension of the sample at least as $`ML^\zeta `$.

Especially in higher dimensions, direct simulations of eq. (2) are not viable, because they violate a crucial no-passing theorem,Middleton  which we respect.

Because of the expected universality and the difficulties of direct simulation, attempts to compute the critical exponents of the quenched EW equation considered discrete systems with non-hamiltonian dynamic rules. These cellular automata models Sneppen ; Buldyrev  are not clearly connected with an equation of motion in the continuum limit Vannimenus  and have often given only rough estimates for $`\zeta `$.

Nevertheless, simulations using cellular automata illustrated that critical elastic manifolds in random media are not generally described by the quenched EW equation, but may fall into two broad classes: The automata Nattermann2 ; Leschhorn  summarized in Ref. Nattermann1, , yield critical exponents which are close to those obtained by analytical work on eq. (2); on the other hand the automata Sneppen ; Buldyrev  summarized in Ref. BuldyrevII,  give very different values for the exponents, incompatible with the results obtained by renormalization group results, but closer to experiment.Buldyrev ; BuldyrevIII ; Krim  There has been much confusion about the actual values of these exponents for $`d>1`$. It is not clear whether the upper critical dimension of this class is also $`d_{\text{uc}}=4`$.

In our controlled Hamiltonian approach, we recover this very rich behavior as a dependence of the roughness exponent $`\zeta `$ on the functional form of the elastic energy: We investigate elastic energies of the form $`E_{\text{el}}|\mathrm{\Delta }|^\gamma `$, and refer to them as $`\mathrm{\Delta }^\gamma `$–models. The $`\mathrm{\Delta }^2`$–model studies directly the eq. (2). Models with $`\gamma 2`$ differ in their critical behavior from the harmonic model. A special role is played by the $`\mathrm{\Delta }^4`$–model, which represents the first non-harmonic corrections of a general elastic energy $`E_{\text{el}}a\mathrm{\Delta }^2/2+b\mathrm{\Delta }^4/12+\mathrm{}`$

We first discuss this important issue for an elastic string ($`d=1`$) in a $`2`$–dimensional medium, which has been abundantly studied in the past. As shown in Fig. 2, we find for the $`\mathrm{\Delta }^2`$–model a value of the roughness exponent $`\zeta _{\mathrm{\Delta }^2}=1.26\pm 0.01`$, which compares well with the two-loop calculation of Chauve et al., while clearly excluding the one-loop result $`\zeta _{\mathrm{\Delta }^2}^{\text{oneloop}}=1`$. A value in excess of one for the $`(1+1)`$–dimensional roughness has been found in direct numerical integrations of the equations of motion,Jensen  in Monte Carlo simulations,RossoKrauth1 ; RossoKrauth2  and in the cellular automata of Ref. Nattermann1, .

Figure 2 also shows our results for the $`\mathrm{\Delta }^4`$–model. We find f-n-inset  $`\zeta _{\mathrm{\Delta }^4}=0.635\pm 0.005`$. This exponent coincides with the value found in Ref. BuldyrevII, . In lattice models, we already showed that the roughness exponent does not change for the $`\mathrm{\Delta }^6`$-model and even in the presence of a metric constraint (bounded $`|\mathrm{\Delta }|`$). In this sense, the $`\mathrm{\Delta }^4`$–model is only one representative of systems with a stronger than harmonic elastic energy. The values obtained from Fig. 2 are much more precise than the earlier ones.

In this problem, the study of rotational invariance of the equations of motion has occupied a crucial role. In fact, eq. (2), and therefore the equation of motion of our $`\mathrm{\Delta }^2`$–model, is invariant under a tilt of the manifold

$$h(𝐱,t)h(𝐱,t)+m_ix_i.$$
(5)
It was observed Amaral  that the automata of Ref. BuldyrevII,  present strong dependence of the velocity and of the value of $`f_c`$ on the parameters $`m_i`$. In our $`\mathrm{\Delta }^4`$–model, rotational invariance is manifestly broken as an elastic energy $`b\mathrm{\Delta }^4/12`$ generates a piece

$$b^2h(h)^2$$
(6)
in the equation of motion. From dimensional analysis,Kardar ; KardarAniso  among the orientation-dependent differential operators, the KPZ term $`(h)^2`$ is more relevant than the term eq. (6). The non-linear KPZ term appears in domain-growth models without disorder, where it is coupled to the velocity.KPZ  Adding it to eq. (2) constitutes the quenched KPZ equation. Its critical exponents have not been computed within the functional renormalization group, and it is not even clear how this term is generated in the absence of a finite velocity, at depinning. Narayan

Kardar et al. Kardar ; KardarAniso  have suggested the force anisotropy due to the disorder as a possible generating mechanism, but none of the automata models present such an anisotropy. However, in their construction rules, a metric constraint is hidden. This implements a strong elastic potential, and naturally generates terms like eq. (6). In our opinion, the main unresolved theoretical issue is to understand how this term reduces to the KPZ term under coarse-graining.

As discussed in Ref. RossoKrauth2, , there is an easy way to understand why the ($`1+1`$)–dimensional $`\mathrm{\Delta }^2`$–model is unstable to higher-order corrections: A roughness $`\zeta >1`$ impliesTangLeschhornDongcomm  that the mean local elongation $`|\mathrm{\Delta }|`$ grows with the system size at least as $`L^{\zeta 1}`$, and thus diverges in the thermodynamic limit $`L,M\mathrm{}`$. Higher-order terms in the elastic energy are thus trivially relevant in one dimension.RossoKrauth2

In $`2+1`$ dimensions, we have studied samples up to size $`L^2\times M=512^2\times 512`$. For the $`\mathrm{\Delta }^2`$–model (shown in Fig. 3), we find excellent scaling for a roughness $`\zeta _{\mathrm{\Delta }^2}=0.753\pm 0.002`$, which again falls between the results of the one-loop and two-loop renormalization group calculations. Our result improves by an order of magnitude the precision of the previous estimate.Nattermann1  Note that our algorithm allows us to know, without invoking universality arguments, that the exponent $`\zeta _{\mathrm{\Delta }^2}0.753`$ is that of the two-dimensional quenched EW model.

In Fig. 3, we also show our results for the $`\mathrm{\Delta }^4`$–model. We find a different value for $`\zeta `$, namely $`\zeta _{\mathrm{\Delta }^4}=0.45\pm 0.01`$. This last value is significantly smaller than $`\zeta 0.48`$ found in Ref. BuldyrevII,  using a cellular automaton.

The two-dimensional $`\mathrm{\Delta }^2`$–model is not only unstable with respect to stronger elastic potentials, but also to any weaker ones. We have studied the $`\mathrm{\Delta }^\gamma `$-model, for $`\gamma `$ slightly below $`2`$: Already for $`\gamma =1.8`$ (cf Fig. 3), we find that the roughness exponent changes drastically, as we obtain clear indications that $`\zeta _{\mathrm{\Delta }^{1.8}}>1`$. This implies that the surface breaks.

We also studied elastic energies $`\sqrt{1+\mathrm{\Delta }^2}`$ in the ($`2+1`$)–dimensional problem. In one dimension, this functional form corresponds to the length of the string, and the $`\mathrm{\Delta }^2`$–model is often understood as the first term from an expansion of the length in powers of $`\mathrm{\Delta }`$. The square root is softer than its first expansion coefficient; we again find values of $`\zeta _{\sqrt{1+\mathrm{\Delta }^2}}`$ in excess of one for large system sizes.

As in $`1+1`$ dimensions, we also find for the ($`2+1`$)–dimensional $`\mathrm{\Delta }^2`$–model indications of a (much slower) divergence of the mean local extension, as shown in Fig. 3. For the $`\mathrm{\Delta }^2`$–model, it seems to diverge logarithmically, whereas for the $`\mathrm{\Delta }^4`$–model it saturates already for very small systems. Again it appears that only the $`\mathrm{\Delta }^4`$–model has a proper thermodynamic limit for the critical manifold.

We also performed extensive computations in $`3+1`$ and in $`4+1`$ dimensions, where we have reached sample sizes of $`L^3\times M=64^3\times 64`$ and $`L^4\times M=32^4\times 32`$, respectively. For the harmonic model in $`3+1`$ dimensions, we obtain $`\zeta _{\mathrm{\Delta }^2}=0.355\pm 0.01`$ again in good agreement with the numerical results of Ref. Leschhorn, , while our result $`\zeta _{\mathrm{\Delta }^2}=0.25\pm 0.02`$ for the $`\mathrm{\Delta }^4`$–model is in contradiction to Ref. BuldyrevII,  which suggested an exponent $`\zeta 0.38`$ larger than for the harmonic case. We checked very carefully that our estimate for $`\zeta `$ is independent of the parameters $`a`$ and $`b`$.

At last, our algorithm is able to investigate $`4+1`$ dimensions.RossoKrauthVannimenus  We have computed manifolds of sizes $`L=8,16,32(L=M)`$. Even though we have only three points, we find for the $`\mathrm{\Delta }^2`$–model an inconsistent fit to a functional form $`W^2(L)L^{2\zeta }`$ for sizes $`L`$ where it was already good in $`3+1`$ dimensions. At the upper critical dimension, one would expect logarithmic behavior for $`W^2(L)`$. The most natural explanation of our findings in $`4+1`$ dimensions is that the $`\mathrm{\Delta }^2`$–model’s upper critical dimension is $`d_{\text{uc}}=4`$. This will still have to be confirmed for larger sample sizes.

For the $`\mathrm{\Delta }^4`$–model, we also do not find a good fit for $`\zeta `$, but the deviations from the functional form $`W^2(L)L^{2\zeta }`$ are less striking than for the $`\mathrm{\Delta }^2`$–model. It is unclear to us whether the roughness exponent of the $`\mathrm{\Delta }^4`$–model *can* be higher than for the $`\mathrm{\Delta }^2`$–model.

In conclusion, we have directly computed roughness exponents for $`d`$–dimensional manifolds in a ($`d+1`$)–dimensional embedding space for $`d=1,\mathrm{},4`$. Some of our findings nicely fit into the existing theoretical framework. For example, the upper critical dimension $`d_{\text{uc}}=4`$, conjectured on the basis of dimensional power counting, is consistent with our data for the $`\mathrm{\Delta }^2`$–model.

In the harmonic case, agreement with the two-loop calculation is fair, and the exactness of the one-loop result can be excluded (cf inset of Fig. 4). Our numerical work relies on no additional hypothesis, as we compute well-defined critical manifolds in finite systems. Even for the harmonic case, though, we would like to understand the divergence of the local scale, for which we find evidence even in $`d=2`$: slightly sub-harmonic potentials lead to a break-up of the critical surface.

For the first time we have presented a Hamiltonian model which reproduces the results of the quenched KPZ class. The connection between our $`\mathrm{\Delta }^\gamma `$–model and the equations of motion in the continuum limit is transparent. The analysis of the $`\mathrm{\Delta }^4`$–model gives a very precise estimate for the value of $`\zeta `$ as a function of dimension $`d`$. It would be very interesting to understand how the local differential operators, as $`^2h(h)^2`$, which correspond to elastic Hamiltonians, renormalize into more relevant terms, as the KPZ term.

Acknowledgments: We thank P. Le Doussal for very helpful discussions all along this work, and also acknowledge stimulating discussions with D. S. Fisher and K. J. Wiese. Our computations were performed on clusters of workstations at Ecole Normale Supérieure, the Institut für Theoretische Physik of the Universität Magdeburg (Germany) and at the Paderborn Center for Parallel Computing (Germany). AKH was financially supported by the DFG (Deutsche Forschungsgemeinschaft) under grant Ha 3169/1-1.


text:
# Quantum comparison machines with one-sided error

## Abstract

It is always possible to decide, with one-sided error, whether two quantum states are the same under a specific unitary transformation. However we show here that it is impossible to do so if the transformation is anti-linear and non-singular. This result implies that unitary and anti-unitary operations exist on an unequal footing in quantum information theory.

The formalism of quantum mechanics has given rise to many theorems concerning ‘impossible machines’. For example, many quantum machines such as quantum copiers  and quantum erasers  are impossible. However these operations can be performed approximately . Furthermore, Galvão and Hardy have shown that approximate cloning can actually enhance the performance of some computations . It is therefore very important to understand the limitations of operations which are ‘approximately impossible’.

Here we consider quantum decision machines with a one-sided error probability. These machines output the answer YES or NO, however the probability of error is non-zero for one of the two possible outputs. Although this concept sounds quite theoretical, practical examples of such one-sided machines are not unfamiliar to physicists. In particular the interaction-free measurement scheme introduced by Elitzur and Vaidman in Ref.  is inherently a decision machine with one-sided error. In this article, we study the fundamental limit of a subclass of quantum decision machines, which we call comparison machines. We show that although it is always possible to decide—with one-sided error—whether two states are the same under a specific unitary transformation, it is impossible to do so if the transformation is anti-linear and non-singular. The present work therefore demonstrates the unparallel roles of unitary and anti-unitary transformations in quantum information theory.

Without loss of generality, we assume that the answer YES or NO is given by a measurement on the first qubit of the output state, and that $`|0`$ corresponds to YES while $`|1`$ corresponds to NO. An immediate example of a decision machine is the SWAP-test discussed in Ref. . Given two physical systems in unknown states $`|\varphi `$ and $`|\psi `$, it is impossible to tell whether they are identical or not—however Buhrman et al have shown in Ref.  that it is indeed possible to do so in the presence of some one-sided error. Explicitly, if $`|\varphi =|\psi `$, then the machine will always say YES. However if $`|\varphi |\psi |=\delta <1`$, then the machine says YES with probability $`(1+\delta ^2)/2`$. This SWAP-test was also utilized by Gottesman and Chuang to construct their digital signature scheme .

A further motivation for studying such a comparison machine, is that the output is a physical observable. For instance, we may use the SWAP-test to determine the payoff to Clare (the cloner) in Werner’s cloning game  instead of the fidelity. For example, we may assign a payoff of 1 to Clare if her clones pass the test, and $`1`$ otherwise. As the SWAP-test depends only on the fidelity of the two states, in terms of expected payoff, this game and Werner’s original one are equivalent. However the game is now physically implementable, whereas Werner’s orignal game was not because there is no physical operation to compute the fidelity of two states. Furthermore, with the SWAP-test adopted as the gauge for payoff, the game formalism and the Quantum Minimax Theorem developed in Ref.  can be applied to show that Werner’s cloning operation is indeed the strategy at the Nash equilibrium .

Given two states $`|\varphi `$ and $`|\psi `$, a more general version of the SWAP-test is to test, with some one-sided error probability, whether $`|\psi =K|\varphi `$ for some arbitrary operation $`K`$. If $`K`$ is unitary, this is made trivially possible using the SWAP-test on $`K|\varphi `$ and $`|\psi `$. On the other hand, we now show that such a device is not possible if $`K`$ is anti-linear and non-singular. Given an arbitrary two-level state: $`a|0+b|1`$, the state orthogonal to it is $`\overline{a}|1\overline{b}|0`$, so an anti-linear map is needed to transform an arbitrary two-level state to its orthogonal state. Therefore, our result immediately rules out the possibility of deciding whether two arbitrary states are orthogonal with one-sided error probability. Our result suggests that it is much more stringent to test the orthogonality condition than to test the identity condition. A complementary idea was discussed in Ref.  where Gisin and Popescu showed that antiparallel pairs of states encode more information than parallel pairs of states. Our result also contrasts with the many proofs of impossibility covering both unitary and anti-unitary maps , and hence demonstrates concretely the unparallel roles of unitary and anti-unitary transformations in quantum information theory.

We now give the proof of the statement, noting that the proof only relies on the linearity assumption of quantum mechanics. Given an arbitrary state $`|\varphi `$ of particle 1 and another arbitrary state $`|\psi `$ of particle 2, let $`M`$ be a machine capable of deciding whether $`|\psi =K|\varphi `$ with some one-sided error probability. We now denote $`K|0`$ by $`|\varphi _0`$ and $`K|1`$ by $`|\varphi _1`$. We first suppose that the machine says NO ($`|1`$) with certainty, but it errs when the answer is YES ($`|0`$):

$`M:|0_1|\varphi _0_2|ancillar`$ $``$ $`a_{00}|0|A_{00}+b_{00}|1|B_{00}`$
$`|1_1|\varphi _1_2|ancillar`$ $``$ $`a_{11}|0|A_{11}+b_{11}|1|B_{11}`$
$`|0_1|\varphi _1_2|ancillar`$ $``$ $`|1|B_{01}`$
$`|1_1|\varphi _0_2|ancillar`$ $``$ $`|1|B_{10}.`$
However, since $`K`$ is non-singular, $`|\varphi _1K(\frac{1}{\sqrt{2}}(|0+|1))`$, hence $`M(\frac{1}{\sqrt{2}}(|0_1+|1_1)|\varphi _1_2)=|1|`$ for some state $`|`$. This implies that $`a_{11}=0`$. Similarly one can show that $`a_{00}=0`$, hence the machine is trivial.

We next suppose that the machine says YES with certainty, but errs when it says NO. In this case

$`M:|0_1|\varphi _0_2|ancillar`$ $``$ $`|0|A_{00}`$
$`|1_1|\varphi _1_2|ancillar`$ $``$ $`|0|A_{11}`$
$`|0_1|\varphi _1_2|ancillar`$ $``$ $`a_{01}|0|A_{01}+b_{01}|1|B_{01}`$
$`|1_1|\varphi _0_2|ancillar`$ $``$ $`a_{10}|0|A_{10}+b_{10}|1|B_{10}`$
Using linearity, $`M(\frac{1}{2}(|0_1+|1_1)(|\varphi _0_2+|\varphi _1_2))=\frac{1}{2}|0(|A_{00}+a_{01}|A_{01}+a_{10}|A_{10}+|A_{11})+\frac{1}{2}|1(b_{01}|B_{01}+b_{10}|B_{10})`$. However since $`K(|0+|1)=(|\varphi _0+|\varphi _1)`$, this implies $`b_{01}|B_{01}+b_{10}|B_{10}=0`$. On the other hand, $`M(\frac{1}{2}(|0_1+i|1_1)(|\varphi _0_2i|\varphi _1_2))=\frac{1}{2}|0(|A_{00}+ia_{10}|A_{10}ia_{01}|A_{01}+|A_{11})+\frac{i}{2}|1(b_{10}|B_{10}b_{01}|B_{01})`$. Again since $`K(|0+i|1)=(|\varphi _0i|\varphi _1)`$ by anti-linearity, we have $`b_{01}|B_{01}b_{10}|B_{10}=0`$. Therefore, the machine always outputs YES. This completes the proof of impossibility.

In conclusion, we have ruled out any physical devices that would decide whether two states are the same under a non-singular anti-linear transformation. As a result, unitary and anti-unitary operations are placed on an unequal footing. This result suggests that further investigation is needed to clarify the boundary of feasible operations in quantum information theory.

CFL thanks NSERC (Canada), ORS (U.K.) and the Clarendon Fund (Oxford) for financial support.


text:
# Rarefaction Shocks, Shock Errors and Low Order of Accuracy in ZEUS

## 1 Introduction

ZEUS is a freely available MHD code that is widely used by the Astrophysical community. Although Stone & Norman (1992a,b) give results for the Sod problem (Sod 1978) and its MHD equivalent, the Brio and Wu problem (Brio & Wu 1988), ZEUS does not appear to have been tested on a wide range of Riemann problems such as those described in e.g. Dai & Woodward (1994), Ryu & Jones (1995), Falle, Komissarov & Joarder (1998)and Balsara (1998).

Since ZEUS is neither upwind for all characteristic fields nor conservative, we might expect it to perform significantly less well than upwind conservative codes (e.g. Brio & Wu 1988; Dai & Woodward 1994; Ryu & Jones 1995; Falle, Komissarov & Joarder 1998; Balsara 1998; Powell et al. 1999). As we shall see, this is indeed true in the sense that there are a number of simple problems for which the ZEUS solution contains significant errors that are absent in solutions calculated with an upwind conservative scheme.

## 2 Rarefaction Shocks

Figures 1 and 2 show that ZEUS generates rarefaction shocks for both pure gas rarefactions and fast magnetosonic rarefactions, whereas the upwind conservative scheme described in Falle, Komissarov & Joarder (1998) gives quite satisfactory results. In both cases, the ZEUS solutions are sensitive to the inertial frame and the rarefaction shocks can be removed by a Galilean transformation that increases the x velocity sufficiently.

These rarefaction shocks are steady structures whose width does not increase with time. Since the effect of the nonlinear terms is to spread such structures, it is clear that the truncation errors in ZEUS must be anti-diffusive in these cases. The most obvious explanation for this is that ZEUS is second order in space, but first order in time since this can lead to an anti-diffusive term in the truncation error. For example, the upwind scheme can be made first order in time and second order in space by omitting the preliminary first order step and in that case it can be shown to be anti-diffusive and also produces rarefaction shocks.

Although ZEUS is second order in space and time for linear advection, the use of a partially updated velocity in the advection step means that it is only first order in time if the velocity is not constant. Further evidence that this is the cause of the problem is provided by the sensitivity of the rarefaction shocks to the Galilean frame and the fact that they disappear when the Courant number is reduced from $`0.5`$ to $`0.1`$, whereas they become much worse if the Courant number is increased above $`0.5`$.

ZEUS has a facility for adding a linear artificial viscosity whose magnitude is determined by the parameter $`qlin`$. The addition of such a viscosity removes the anti-diffusive terms by reducing the scheme to first order in space for everything except linear advection. For the gas rarefaction, $`qlin=0.25`$ cures the problem and seems to be optimal for a global Courant number of $`0.5`$, but it is too large if the local Courant number associated with the wave is small. Since the linear viscous term must balance an anti-diffusive term that scales like the timestep, it would be better if the viscous term that is implemented in ZEUS were multiplied by the local Courant number associated with the wave that is causing the problem. Since rarefaction shocks only arise for rarefactions associated with the sound wave with the largest speed relative to the grid, it is the smallest local Courant number that is appropriate.

In MHD the situation is even worse since, although the rarefaction shocks in the fast rarefaction can be removed by setting $`qlin=1`$, this makes the scheme very diffusive for other waves. Furthermore, the required value of $`qlin`$ depends on the particular problem. It might be possible to avoid such a large value of $`qlin`$ by adding an appropriate artificial resistivity, but the code has no facility for this.

Figures 3 and 4 show that, even for an initially smooth rarefaction wave, ZEUS is significantly less accurate than an upwind scheme. The results are for $`qlin=0.25`$, but Figure 4 shows that ZEUS is still first order even without this. In contrast, it is evident from Figure 4 that the rate of convergence of the upwind scheme is second order. Note that the upwind scheme also has an artificial viscosity as described in Falle, Komissarov & Joarder (1998), but since this is applied in the Riemann solver, it does not reduce the order in smooth regions. Incidentally, ZEUS performs even worse if one does not take the staggered grid into account in setting up the the initial solution. Furthermore, for both codes, point samples were used to project the exact solution onto the grid, which is reasonable for ZEUS, but is somewhat unfair to a conservative scheme.

The upwind scheme produces reasonable results at the lowest resolution, even though this corresponds to only $`2`$ cells in the rarefaction at the initial time, whereas ZEUS needs at least $`4`$ cells for the same accuracy. For a three dimensional calculation, this would require $`24`$ times the computing time and $`8`$ times the memory since ZEUS is about $`2/3`$ times the speed of the upwind scheme. The disparity in efficiency is actually greater than this because for both codes the Courant number was set to the ZEUS default value of $`0.5`$ for all cases described in this paper. The upwind code can run at larger Courant numbers than this, whereas even $`0.5`$ can be too large for ZEUS for some Riemann problems. Of course, the slower convergence of ZEUS also means that the situation would be even worse if greater accuracy were required.

## 3 Shock Errors

Since ZEUS is not conservative, we expect it to generate errors at shocks which cannot be reduced by increasing the resolution. As it turns out, these errors are small ($`<5\%`$) for pure gas dynamics and are entirely absent for an isothermal equation of state. However, they can be significant for adiabatic MHD .

Figure 5 shows that, for a nearly perpendicular fast shock, the post-shock gas pressure in the ZEUS solution is too low by a factor of 2. In contrast, the conservative upwind scheme gets the solution exact to rounding. It is true that this is a somewhat extreme case since the plasma $`\beta `$ is negligible upstream of the shock and $`\beta =0.037`$ downstream. However, such low values of $`\beta `$ do occur in dense molecular clouds and protostellar discs (e.g. Crutcher 1999). Furthermore, even though $`\beta `$ is small, such errors in the gas pressure can have a significant effect on the dynamics because the gas pressure provides a force parallel to the field, whereas the Lorentz force does not.

Finally, Figure 6 shows that a relatively small error at a fast shock can be amplified by a slow shock following on behind. In this case the ZEUS solution has an error of $`22\%`$ in the density behind the slow shock travelling to the right. This is not caused by small $`\beta `$ since $`\beta =0.16`$ behind the fast shock, $`\beta =6.1`$ behind the slow shock and the error in the gas pressure is much smaller than in the density.

Like Balsara (2001), we find that ZEUS produces large post-shock oscillations for strong MHD shocks, but that these can be reduced substantially by adding the same linear artificial viscosity that removes gas dynamic rarefaction shocks. This is presumably because a quadratic viscosity leads to algebraic decay of these oscillations, whereas a linear viscosity gives exponential decay. The calculation shown in Figure 6 used this value of the linear artificial viscosity and it can be seen that the amplitude of the post-shock oscillations is quite small.

The calculations presented are all coplanar ($`v_z=B_z=0`$), but we have also looked at some non-coplanar problems in order to see whether the presence of Alfvén waves causes any additional difficulties for ZEUS. This is not the case, at least for the problems that we have considered.

## 4 Conclusion

It is evident from these results that, ZEUS can be made just about acceptable for pure gas dynamics if the linear artificial viscosity is multiplied by the smallest local Courant number since the shock errors are small in this case. However, it is not satisfactory for adiabatic MHD, at least in its present form. The shock errors do not occur for an isothermal equation of state, but, since the rarefaction shocks do, ZEUS is also not reliable for isothermal MHD. It is possible that the rarefaction shocks in MHD waves can be removed without using an excessive linear artificial viscosity by the addition of an appropriate linear artificial resistivity. The shock errors might also be reduced by advecting the total energy rather than the internal energy. However, even with such improvements, the low order of accuracy makes ZEUS very inefficient compared with a modern upwind scheme.

This should not be taken to mean that conservative upwind codes are in any sense perfect. For example, it is necessary to introduce some extra dissipation in the Riemann solver to remove the serious errors discussed by Quirk (1994) and some desirable properties, such as strict conservation, may have to be sacrificed in order to satisfy the constraint $`𝐁=0`$ in multidimensional MHD (see e.g. Powell et al. 1999; Balsara 2001).

These results obviously have implications for the reliability of the numerous calculations in the literature that have used ZEUS. Although these effects are likely to be present in many cases, the associated errors are not necessarily so serious as to completely invalidate the calculations. Whether or not they make any qualitative difference in any particular case can only be decided either by a thorough examination of the results to see whether any of these errors are present, or by repeating the calculations using a modern code.

These calculations were performed with the version of Zeus2d available from the NCSA website, but, since all versions of ZEUS appear to use the same algorithms, the results should not depend on the particular version. It is also worth pointing out that although we used the scheme described by Falle, Komissarov and Joarder (1998), similar results would probably have been obtained with any modern upwind code.

The author would like to thank both the editor and an anonymous referee for a number of helpful comments on the original version.