text:
# BURSTS OF EXTENSIVE AIR SHOWERS: CHAOS VS. STOCHASTICITY

## 1 Introduction

We have already studied the distribution of arrival times of extensive air showers (EAS) registered with the EAS–1000 prototype array both by methods of classical statistics  and by methods of cluster analysis . In particular, we presented 20 EAS cluster events—groups of consecutive showers that were registered in time intervals much shorter than expected ones. These groups represent bursts of the EAS count rate and may be considered as cosmic ray bursts, see . This phenomenon has put forward at least two questions: one on the astrophysical nature of the process, another on its statistical properties. Namely, while the vast majority of sufficiently long samples satisfy the hypothesis for an exponential distribution of time delays between EAS arrival times , the $`\chi ^2`$-test performed for samples taken in the vicinity of some of the EAS clusters made us reject this hypothesis. Thus we have decided to apply methods of nonlinear time series analysis to samples that contain EAS clusters in order to clarify dynamical reasons of this situation. This approach has already proved to be a powerful tool of investigation in different fields of science including different fields of astrophysics in general  and cosmic ray physics in particular, see  and references therein.

Recall that the data set under consideration represents 203 days of regular operation of the array for the period August 30, 1997 to February 1, 1999. The total EAS number in the data set equals 1 668 489. The mean number of charged particles in a shower is of the order of $`1.2\times 10^5`$. This corresponds to the energy of a primary particle $`1`$ PeV/nucleon. The mean interval between consecutive EAS arrival times is equal to 10.5 s. The discreteness in the moments of EAS registration approximately equals 0.055 s (one tic of the PC clock).

In , we have already discussed briefly the results of nonlinear analysis of one of the EAS clusters. In the present paper, we study this cluster in depth and briefly discuss three other clusters that demonstrate nonlinear features.

## 2 Preliminary Facts

In what follows, we study a scalar time series of the form $`x_1,x_2,\mathrm{},x_n`$, where $`x_i=t_it_{i1}`$, $`t_i`$ is the moment of registration of the $`i`$th EAS, $`i=1,\mathrm{},n`$, and $`t_0=0`$ corresponds to midnight. For the purposes of our investigation, it is convenient to use time intervals $`x_i`$ given in seconds.

One can point out two main approaches to nonlinear analysis of the time series: (1) the scaling exponent and fractal length estimates based on the self-similarity properties of experimental signals ; (2) an estimate of the correlation dimension based on the embedding of the original time series in an $`m`$-dimensional phase space  (see, e.g.,  for a comprehensive review). In the current research, we follow the second approach. Within it, the main tool to analyse the dynamics of the time series is the Grassberger–Procaccia method  with the modification suggested by Theiler . Namely, given a sample $`(x_1,x_2,\mathrm{},x_N)`$ extracted from the time series we construct $`m`$-dimensional delay vectors

$$𝐱_i=(x_i,x_{i+\tau },x_{i+2\tau },\mathrm{},x_{i+(m1)\tau }),$$

where $`\tau `$ is an arbitrary but fixed time increment and $`m`$ is an embedding dimension. Then we compute the number $`K(\rho )`$ of vectors with mutual distance less or equal than $`\rho `$ and such that delay vectors $`𝐱_i`$ are shifted by at least $`W`$ indices:

$$K(\rho )=\underset{i=1}{\overset{MW}{}}\underset{j=i+W}{\overset{M}{}}\mathrm{\Theta }(\rho 𝐱_i𝐱_j),$$
(1)
where $`\mathrm{\Theta }`$ is the Heaviside step function, $`M=N(m1)\tau `$, and $`W1`$ is the cut-off parameter (the Theiler window). Finally, we plot $`\mathrm{log}K(\rho )`$ vs. $`\mathrm{log}\rho `$. For small $`\rho `$, the slope of this plot is an estimate of the correlation dimension $`D_2`$:

$$D_2(\rho )=\frac{\mathrm{log}C_2(\rho )}{\mathrm{log}\rho },$$

where $`C_2`$ is the correlation sum:

$$C_2(\rho )=\frac{2K(\rho )}{(MW)(MW+1)}.$$

(Obviously, one may use $`K`$ instead of $`C_2`$ in the expression for the correlation dimension.) A plateau observed in the $`D_2(\rho )`$-plot for small $`\rho `$ or, equivalently, a so called scaling region in the $`\mathrm{log}K`$ vs. $`\mathrm{log}\rho `$ plot are regarded as signs of chaotic dynamics in the corresponding time series. A value of $`D_2`$ at the plateau is taken as an estimate of the correlation dimension of the attractor underlying the data. This quantity also gives a (lower) estimate for the number of degrees of freedom in the process under consideration.

To compute $`K`$, we normally divide each unit interval in $`\mathrm{lg}\rho `$ into 50 subintervals of equal length. We have found that though this number of subintervals seems to be small, it adequately reflects the qualitative structure of $`D_2`$ while a much bigger number of subintervals ($`200`$) leads to considerable fluctuations in $`D_2`$; these fluctuations hide the structure of $`D_2`$, especially for $`m`$ large enough. A smaller number of intervals gives too rough structure of $`D_2`$. The derivative is calculated via a standard three-points algorithm. No smoothing or fitting procedures are used.

Notice that Eq. (1) contains four free parameters: $`N`$, $`\tau `$, $`W`$, and $`m`$. At the preliminary stage of the investigation, we “scanned” the experimental data set having split it into adjacent samples with $`N=128`$, 256, 512, and 1024. At this step, we used $`\tau =1`$, $`W=1`$, and odd values of $`m`$ in the range from 5 to 13. The value of $`N`$ was chosen to be a power of 2 because this allows one to use fast algorithms of calculating the Fourier transform that is the main part of traditional Fourier analysis. After we had found a number of samples that demonstrated some kind of plateau in the plot of $`D_2(\rho )`$, we studied the corresponding data with different $`N`$ in the above range and different values of $`\tau `$ in the range 1–10. (Evidently, in our case $`\tau `$ may take only integer values.) For $`N`$ large enough, we have performed calculations for $`m`$ up to 25. In any case, $`m`$ and $`\tau `$ were chosen such that the number of delay vectors was greater than 100. Finally, to avoid autocorrelation in the time series, we employed $`W=1\mathrm{}20`$.

To compute mutual distances between delay vectors (see Eq. (1)), we tried several norms: the maximum norm $`L_{\mathrm{}}`$, the taxicab norm $`L_1`$, and the Euclidian norm $`L_2`$. Recall that for a given $`m`$-dimensional vector $`𝐱`$ these norms are defined in the following way:

$$𝐱_{\mathrm{}}=\underset{1im}{sup}\left|x_i\right|,𝐱_1=\underset{i=1}{\overset{m}{}}\left|x_i\right|,𝐱_2=\sqrt{\underset{i=1}{\overset{m}{}}\left|x_i\right|^2}.$$

Besides these, we have also employed the “dimension scaled” norms $`L_{1C}`$ and $`L_{2C}`$ that are expected to provide more reliable results than their classical counterparts :

$$𝐱_{1C}=\frac{1}{m}𝐱_1,𝐱_{2C}=\frac{1}{\sqrt{m}}𝐱_2.$$

Since the process of calculating $`D_2`$ is very time consuming, we have used only one of these norms, namely the maximum norm, to analyse the complete data set. Four other norms were tested for a number of samples both with and without a plateau in the $`D_2(\rho )`$-plot obtained with $`L_{\mathrm{}}`$.

There are a number of tools that can help one to verify the results of calculating the correlation dimension $`D_2`$, see, e.g., . Among them, we have chosen the Theiler–Takens “maximum likelihood” estimator of the correlation dimension :

$$t_{\mathrm{TT}}(\rho )=C_2(\rho )\left[_0^\rho \frac{C_2(\rho ^{})}{\rho ^{}}𝑑\rho ^{}\right]^1.$$
(2)
This quantity not only provides an efficient means to find out an optimal estimate of the correlation dimension out of the correlation integral but also provides an estimate of the statistical error.

Another problem in the case when a plateau in a plot of $`D_2(\rho )`$ is observed is to make an assessment about the nature of the dynamics. The main difficulty is to figure out whether one witnesses chaotic dynamics in a deterministic process or a special class of stochastic processes . One of the main tools to solve this problem is the method of surrogate data , see also  and references therein. The main idea of this approach is to generate a sufficient number of artificial (“surrogate”) samples that have the same statistical distribution (and possibly some other features) as the experimental data and to study the behaviour of the correlation dimension. If surrogate data demonstrate a plateau in plots of $`D_2(\rho )`$ then one concludes that the original data is stochastic; otherwise they are likely to be deterministic and chaotic.

Among other tools one can find the quantity

$$t_{\mathrm{BDS}}(m,\rho )=\frac{C_2(m,\rho )}{C_2(1,\rho )^m}$$
(3)
suggested in  and in the above form—in . As it was shown in , for a sequence of independent random numbers, $`C_2(m,\rho )=C_2(1,\rho )^m`$ holds, where $`m`$ is the embedding dimension, and thus $`t_{\mathrm{BDS}}(m,\rho )1`$. Therefore, if one finds $`t_{\mathrm{BDS}}(m,\rho )1`$ for $`\rho `$ where a plateau in the $`D_2(\rho )`$-plots is observed then one may conclude that the nature of the dynamics is chaotic.

Besides this, we used a function suggested in :

$$\varphi _0(\rho )=\frac{D_2(m_2,\rho )D_2(m_1,\rho )}{m_2m_1},$$
(4)
where $`m_1`$ and $`m_2`$ are different embedding dimensions. For independent random data, $`\varphi _0(\rho )0`$ and $`\varphi _01`$ as $`\rho 0`$. We also employed the “normalized slope” introduced in :

$$\varphi (m,\rho )=\frac{1}{m}D_2(m,\rho ).$$
(5)
If $`\varphi (m,\rho )`$ does not converge to 0 in a wide range of $`\rho `$ as $`m`$ grows but to some value $`0.1`$, then most probably data do not represent a chaotic process but should be treated by statistical techniques .

Recall that before applying nonlinear techniques to the analysis of a time series, it is strongly suggested to check whether the time series is really nonlinear. In particular, one can try a measure for time-reversibility, which is considered to be a good indicator for nonlinearity . For the data sorted in time order,

$$\gamma =\frac{1}{(\sigma ^2)^{\frac{3}{2}}(N1)}\underset{i=2}{\overset{N}{}}\left(\frac{x_ix_{i1}}{s_is_{i1}}\right)^3$$
(6)
is calculated, which is just the mean of the slopes, taken to the third power; here $`s_i`$ are the moments of time and $`\sigma ^2`$ is the variance of the sample. (In our case, $`s_is_{i1}=1`$ since $`x_i`$ represent just numerated time delays.) For a time series generated by a linear process and for the surrogate data, one expects $`\gamma 0`$. In contrast, time series with nonlinearities can be asymmetrical in time and may yield values of $`\gamma 0`$. To check, whether $`\gamma `$ significantly deviates from zero for the studied sample, one should generate a sufficient number of surrogate data. To pay regard to deviations in both directions ($`\gamma >0`$ and $`\gamma <0`$), a two sided test has to be performed .

Finally, a few words are in order about stationarity of the time series under consideration. As is well known, a fundamental assumption underlying almost all existing linear and nonlinear techniques of time series analysis is that the time series is stationary, see, e.g., . As we have already mentioned earlier , the count rate of EAS depends on the value of the atmospheric pressure. For the data obtained with the EAS–1000 prototype array, this dependence can be approximately expressed by a simple formula

$$\mathrm{ln}N_{\mathrm{EAS}}=\beta P+\mathrm{const},$$

where $`N_{\mathrm{EAS}}`$ is the number of EAS registered in a time unit, $`P`$ is the atmospheric pressure, mm Hg, and $`\beta 10^2`$ is the barometric coefficient. This effect makes the time series non-stationary. To provide stationarity (at large time scales) we adjusted time delays between consecutive showers to the atmospheric pressure $`P^{}=742`$ mm Hg which is close to the average pressure for the whole analyzed data set. This adjustment was made by the following formula:

$$x_i=x_i^0\mathrm{exp}[\beta (P^{}P_i)],$$

where $`x_i^0`$ is the experimental time delay, and $`P_i`$ is the atmospheric pressure at the $`i`$th delay. The barometric coefficient was chosen to be $`\beta =1.0810^2`$, exactly as in our previous paper . Due to the adjustment, the mean and the variance of the time series at large time scales were approximately constant.

###### Remark 1

In fact, since we shall discuss only comparatively short samples, the above adjustment to $`P^{}`$ is not important for the present results because atmospheric pressure does not change significantly during periods of time covered by the samples. We perform the adjustment in order to guarantee stationarity of the whole time series.

To compute the correlation sum $`C_2`$ effectively, we have worked out an algorithm based on preliminary sorting of mutual distances between delay vectors. Though quite simple, the algorithm occurred to be up to 10 times faster than a straightforward computation of $`C_2`$. To perform calculations, we employed GNU Octave  running in Mandrake Linux.

## 3 The Main Results

As is well known, one needs sufficiently long samples to perform a successful time series analysis. On the other hand, the bursts of the EAS count rate have comparatively short time range. Thus we did not in fact expect to find signs of chaotic dynamics in the vicinity of EAS clusters. Surprisingly, we have found some.

Having “scanned” the available data set, we found that for the overwhelming majority of data samples no scaling region is observed in $`\mathrm{log}C_2`$ vs. $`\mathrm{log}\rho `$ plots. This is quite natural since arrival times of extensive air showers represent a simple stochastic process such that for sufficiently long samples, the number of EAS registered in a time unit obeys the Poisson distribution in a wide range of time delays, see, e.g., . Still, we have come across a number of samples with a plateau in the $`D_2(\rho )`$-plot in the vicinity of four EAS clusters, namely those registered on May 14, November 11, and December 28, 1998 and on January 8, 1999. Let us begin our discussion with the second of these events, which is the longest one.

The cluster event observed on November 11, 1998 consists of 136 EAS registered within the period from 01:21:17.47 to 01:38:02.27 (Moscow local time, MSK111Recall that MSK is connected to UT as follows: $`T_{\mathrm{MSK}}=T_{\mathrm{UT}}+3+\mathrm{\Delta }`$ (hours), where $`\mathrm{\Delta }=1`$ during Daylight Saving Time, 0 otherwise.) with the atmospheric pressure $`P=752.3`$ mm Hg, see Fig. 1. Within the whole data set obtained on November 11, 1998 the showers that form the cluster event have numbers 435–570. The event is made up of three clusters, which begin at consecutive arrival times (i.e., the first cluster begins at the shower #435, the second one begins at the shower #436, etc.) and end up simultaneously at the shower #570. In our opinion, the appearance of three clusters does not reflect any process of astrophysical nature but is caused by the technique of their selection (see  for the details). Thus we treat the event as a single (outer) cluster. The real duration of this cluster equals 1004.8 s while the adjusted duration equals 898.8 s; the probability of the appearance of such a cluster is of the order of $`2\times 10^7`$.

Let us take a look at a number of samples in the vicinity of the cluster in order to see how the correlation dimension $`D_2`$ changes when the cluster appears. To make this influence more clear we shall not only present samples with a plateau in the plot of $`D_2(\rho )`$ but also a typical sample outside the cluster. Besides nonlinear tools, we shall employ the classical Fourier analysis. Recall that it is suggested for the Fourier analysis to have $`x(1)x(N)`$. All samples discussed below are chosen to satisfy this demand. In particular, we omit the last shower of the cluster and consider a sample that consists of 134 instead of 135 delays. Obviously, this does not influence $`D_2`$ and other quantities significantly.

Figure 2 shows the correlation dimension (the left column) and the power spectrum density (PSD) (the right one) calculated for three samples in the vicinity of the cluster. The top row represents a sample that consists of EAS #435–569 and thus contains the cluster without the last delay. The length of the sample $`N=134`$. The middle row represent a sample that contains the cluster but consists of $`N=256`$ delays (EAS #425–681). The bottom row represents a sample that immediately follows the second one. It also consists of 256 delays (EAS #682–938). The correlation dimension was computed for $`\tau =1`$, $`W=1`$, and $`m=1,\mathrm{},12`$.

What can be seen from this figure? In our opinion, the most important thing is that one can see plateaus in both plots of $`D_2(\rho )`$ that represent samples with the cluster. To the contrary, no plateau is observed for the sample outside the cluster. We stress that in this sense, the sample is typical for the whole data set.

Next, it is remarkable that the influence of the cluster on the behaviour of the correlation is not restricted to the cluster itself, e.g., the sample made of EAS #425–681 begins approximately 2.5 min before the cluster and ends up in 23.5 min after it. Our analysis has revealed that a plateau can be observed even for samples with $`N500`$ lying the in the range of EAS numbers 340–1000. It is also interesting to mention that the plateau is more pronounced if a sample contains more of the after-cluster showers than those arrived before.

The level of the plateau changes for different $`N`$. The correlation dimension fluctuates around $`D_21.77`$ for the sample with $`N=134`$ while $`D_22.45`$ for the sample with $`N=256`$. This may be due to low numerousness of the samples and/or white noise effects.222We thank an anonymous referee for pointing this out to us. Notice that for both samples with a plateau the demand $`m>2D_2+1`$ is satisfied if we assume that $`D_2`$ is given by the level of the corresponding plateau (see, e.g., ). We also remark here that the level of a plateau varies for samples with the same length but with different positions with respect to the cluster. The left panel in Fig. 3 shows the dependence of the mean value of $`D_2`$ at a plateau on the position of a sample. All samples have the same length $`N=134`$ but begin at different EAS from #425 up to #505. EAS #425 arrived at 01:18:42.04 MSK, EAS #505 arrived at 01:30:42.15, thus the left ends of these samples cover 12-minute time interval. It is evident from the plot that the value of the correlation dimension takes different values and monotonically decreases from $`1.80`$ to $`1.08`$. For samples that begin at EAS $`510`$ a plateau becomes shorter and soon disappears.

The right panel in Fig. 3 demonstrates how the value of $`D_2`$ saturates as the embedding dimension $`m`$ grows. The plot is made for the sample that consists of EAS #425–681. It is clearly seen that the saturation takes place at $`m=9`$ but not at $`m=3`$ as one may expect for this value of $`D_2`$ and a purely chaotic time series. As it was shown in , this may be due to the presence of observational noise and small length of the sample. The demand $`N10^{D_2/2}`$  is also obviously hold.

Finally, one can notice that the plateaus are not as “smooth” and long as one may like to see. In this connection, we would like to mention two things: (1) Experimental data are affected by noise that can spoil the picture; we plan to employ special noise filtering techniques in future. (2) We do not use any smoothing techniques while calculating $`D_2`$. They can improve a plot of $`D_2`$ but the situation remains qualitatively the same.

Now let us turn to the right column of Fig. 2. This column presents the results of the Fourier analysis of the three samples. The PSD of the sample without a plateau is a broadband one, which is typical for a random noise signal. The PSDs of two other samples do not also differ significantly from a broadband one. Still it seems worth mentioning that for both samples with a plateau one can see the highest peaks located at the left end of the spectrum. This is true for almost all samples with a plateau in $`D_2`$-plots we have found.

Next, let us employ the tools discussed in Section 2 in order to obtain a deeper insight in the observed phenomenon. For brevity, let us call the samples with a plateau in $`D_2`$-plots ‘possibly chaotic’ (PC). Our aim is to figure out whether they are indeed chaotic or represent a special stochastic process.

As one can see from Fig. 4, the Theiler–Takens estimator (2) does not have a clear plateau for either samples but the behaviour of $`t_{\mathrm{TT}}`$ for the PC samples noticeably differs from that for the third one. Namely, for the PC samples the $`t_{\mathrm{TT}}`$-curves saturate for the corresponding intervals of $`\mathrm{lg}\rho `$ for $`m>8`$ while for the third sample the curves stay separated.333We would like to mention that one can find samples such that sufficiently clear plateaus can be observed for both $`D_2`$\- and $`t_{\mathrm{TT}}`$-plots, see the first version of this paper for the details .

Next, the right column of Fig. 4 depicts the behaviour of $`t_{\mathrm{BDS}}`$, see Eq. (3). One can see that for both PC samples and $`m`$ big enough, $`t_{\mathrm{BDS}}`$ noticeably differs from 1 at the intervals where plateaus of the correlation dimension are observed \[$`\mathrm{lg}\rho (0.94,1.12)`$ for the first sample and $`\mathrm{lg}\rho (0.94,1.08)`$ for the second one\]. Conversely, $`t_{\mathrm{BDS}}1`$ for $`\mathrm{lg}\rho 0.94`$ for the third sample. As it follows from the properties of $`t_{\mathrm{BDS}}`$ discussed above, this kind of behaviour witnesses in favour of a chaotic nature of the PC samples and a stochastic nature of the sample without a plateau in the $`D_2`$-plot.

Fig. 5 depicts the normalized slope $`\varphi `$ (see Eq. (5)) and the function $`\varphi _0`$ (Eq. (4)) for the same three samples. One can see that the behaviour of these functions for the PC samples noticeably differs from that for the third one. Still, the situation is not unequivocal. Namely, for both PC samples, $`\varphi _0`$ fluctuates around zero at the corresponding intervals of $`\mathrm{lg}\rho `$ thus giving an argument in favour of a chaotic nature of these samples. On the other hand, the limiting values of the normalized slope $`\varphi `$ (the bottom curves) lie above 0.1 at the same intervals of $`\mathrm{lg}\rho `$ thus suggesting a stochastic nature of the samples. This controversy can probably be resolved if we recall that $`\varphi `$ converges slowly, see . Really, the level of the “hollows” in the plots of $`\varphi `$ decreases if we increase $`m`$ reaching the values $`0.1`$ and $`0.15`$ for the first and the second PC samples respectively and $`m=15`$. For $`m>15`$ the number $`K`$ of close delay vectors at the corresponding intervals becomes $`100`$ thus not providing a value statistically necessary for accurate computation of the correlation dimension.

Next, let us say a few words about the results of surrogate data tests. For the first PC sample, we generated 99 “shuffled” surrogates obtained by a random permutation of time intervals $`x_i`$ in a given sample. This number of surrogates corresponds to a 98% level of significance (L.S.) of the statistical test . For the second sample, we made 39 “shuffled” surrogates (95% L.S.) and 39 surrogates based on the amplitude adjusted Fourier transform method proposed in . To obtain these surrogates, we employed the TISEAN package . None of these “Fourier-based” surrogates demonstrated a plateau in $`D_2`$-plots thus revealing that the consequence of time delays that constitute the original data is crucial for the appearance of plateaus. Therefore, the results of this test suggest that the PC samples represent a chaotic (and thus nonlinear) process.

On the other hand, the test for time-reversibility based on $`\gamma `$-statistic (6) did not detect nonlinearity in the PC samples. Namely, for the first PC sample $`\gamma =0.74`$ while $`\gamma (1.38,1.28)`$ for the surrogates. For the second PC sample, $`\gamma =0.22`$ while surrogates gave $`\gamma (0.44,0.45)`$. (To compare, $`\gamma =0.38`$ for the third sample.) These results imply that the hypothesis for time-reversibility of the PC samples cannot be rejected. Therefore, these samples may be linear (or obtained by a transformation of a linear process) thus excluding the possibility of chaotic dynamics.

Thus, it is interesting to figure out what kind of stochastic processes can describe the distribution of time delays in the PC samples. To do this, we performed the $`\chi ^2`$-test to verify the hypothesis that time delays between consecutive arrival times have an exponential distribution. We used different time bins in the range from 1 to 10 s with a step equal to 0.5 s providing that each bin taken into account contains more than 10 events. We have found that for both PC samples, this hypothesis should be rejected with at least 90% L.S. while for the third sample the hypothesis may be accepted with 85% L.S.

In addition, we remark that for the above PC samples a plateau was only found for $`\tau =1`$ and the maximum norm. Still there are PC samples in the vicinity of the cluster that demonstrate a plateau in the $`D_2`$-plots for $`\tau =2`$ and both “dimension scaled” norms, see  for the details. For all PC samples, a plateau is observed in a wide range of the Theiler window $`W`$ (up to 20).

Now let us briefly discuss three other EAS clusters that produce signs of chaotic behaviour in the corresponding time series. Table 1 contains certain information about these clusters as well as the cluster discussed above. Notice that while the cluster registered on January 8, 1999 is sufficiently long ($`N_{\mathrm{EAS}}=134`$) two other clusters are short and this makes an appearance of signs of chaos in much longer samples very surprising. It is also worth mentioning that the fourth cluster (December 28, 1998) was registered within the period of observation of the GRB No. 7285 in the BATSE Catalogs .

Table 2 presents the results of calculation of the correlation dimension for a number of samples in the vicinity of the clusters as well as the results of a number of tests discussed above. For the last three samples, the given values of $`\overline{D}_2`$ were computed for $`m=15`$, 17, and 14 respectively. Notice that the PC sample for January 8, 1999 is shifted with respect to the cluster. The cluster by itself does not produce signs of chaotic dynamics. For this event, a plateau can be observed not only for $`\tau =1`$ but also for $`\tau =2`$ and 3. For the cluster registered on May 14, 1998, a plateau can also be observed for $`\tau =2`$. For the last cluster, a plateau was observed only for $`\tau =1`$. To the contrary to most other PC samples, a plateau exists not only for the maximum norm $`L_{\mathrm{}}`$ but also for the norms $`L_2`$ and $`L_{2C}`$. For all PC samples, a plateau was found in a wide range of values of the Theiler window $`W`$. This means that autocorrelation, which can lead to spurious results, is avoided.

As one can see from Table 2, the results of the tests applied to these three events are mostly similar to those for the first cluster. Namely, none of the surrogate samples (39 “shuffled” and 39 Fourier-based surrogates) demonstrates a plateau in the $`D_2(\rho )`$-plot thus implying a chaotic nature of the original samples. On the other hand, the test for time-reversibility usually does not detect nonlinearity in the PC samples. The only exception is the last sample in Table 2. For this sample, the $`\gamma `$-test applied to shuffled surrogates gives the usual result but being applied to Fourier-based surrogates it rejects the hypothesis for time-reversibility with 95% L.S. It is worth mentioning that this result has been confirmed when we employed another technique of making Fourier-based surrogates also implemented in TISEAN. Next, for all PC samples the Theiler–Takens estimator $`t_{\mathrm{TT}}`$ demonstrates either a plateau or an interval where the curves saturate as the embedding dimension $`m`$ grows. Finally, it is interesting to note that for three PC samples we have found a way of grouping time intervals that allows one to accept a hypothesis for an exponential distribution of time delays within the corresponding PC sample. Still we must remark that for the majority of groupings this hypothesis should be rejected with at least 90% L.S.

Thus we come to a conclusion that the observed dynamics of time intervals between consecutive EAS in the vicinity of these four clusters may indeed be chaotic though this cannot be assessed unequivocally at the moment. Further investigation is necessary.

## 4 Discussion

To get a deeper insight into physics of EAS clusters, let us consider the cluster registered on Nov. 11, 1998 from another point of view. A natural question that arises when one finds a burst of EAS count rate is whether primary particles that led to the appearance of the burst had had a common compact source. Fig. 6 shows arrival directions of EAS that belong to the cluster and have zenith angle $`\theta <40^{}`$. Obviously, one cannot give a positive answer to the above question. Still, this seems quite natural since the energy of primary particles in our case is too low to “remember” the source. Besides this, the majority of EAS in the cluster represent a normal count rate. Namely, since we register 5.7 EAS/min in average and the duration of the cluster adjusted to $`P^{}=742`$ mm Hg equals approximately 15 min then we may expect to register 85–86 EAS within this period. Naturally, arrival directions of these EAS should be distributed more or less isotropically. Thus, one may expect to observe a kind of a cluster in arrival directions for only remaining 50 EAS. Still, this does not happen. One can see doublets and triplets of showers with almost coincident arrival directions and a kind of grouping along the Galactic plane but no big groups of EAS are observed.

Another question is whether clusters consist of EAS that are more “powerful” than the others. Surprisingly enough, but the answer is negative. As we have already mentioned in Section 1, the mean value of electrons $`\overline{N}_e^{\mathrm{tot}}`$ in EAS in our data set is of the order of $`1.2\times 10^5`$ particles. To compare, for the above cluster $`\overline{N}_e9.7\times 10^4`$, and only 14 EAS have $`N_e>\overline{N}_e`$, see Fig. 6. Two of these EAS have $`N_e10^6`$. On the other hand, at least 41 EAS have $`N_e<1/2\overline{N}_e`$. (We remark that the geometry of the EAS–1000 prototype array does not allow one to obtain parameters of all registered EAS.)

A similar situation is observed for other EAS clusters, both PC and “ordinary” ones. Thus, we come to the following observations: (i) EAS in a cluster do not have a joint source (arrival direction); (ii) the majority of EAS in a cluster have $`N_e<\overline{N}_e^{\mathrm{tot}}`$. This allows us to suggest a simple conceptual model that is intended to explain qualitatively how EAS clusters might be produced. To do this, we need only to assume that a considerable part of EAS that constitute clusters are generated by charged particles.

As is well known, the interstellar space is filled with strong and highly inhomogeneous magnetic fields. Still, average magnetic fields in the heliosphere are not strong enough to influence the dynamics of particles with $`E10`$ TeV/nucleon (see, e.g., ). On the other hand, a series of long-term investigations based on data obtained with Voyager 1, 2, and other space crafts have demonstrated that the large-scale magnetic field strength fluctuations frequently have large amplitudes and are intermittent, and that regions of relatively intense magnetic fields can have a radial extent of more than 10 AU, see  and references therein.

Basing on these results, let us consider the following situation. Suppose there is an extended and sufficiently “thick” layer of accelerated particles (LAP) (possibly called a “wave”) that moves through space towards Earth, see Fig. 7. Normally, if it just passes through, an EAS array registers only a few particles originated from the LAP. Now suppose that the LAP meets a region of an extended strong and inhomogeneous magnetic field that works as a lens. If it happens that this “magnetic lens” declines particles in way that they get focused in a “proper” direction then an array may register an excess of EAS over the normal count rate, i.e., an EAS cluster. Notice that from an observer’s point of view showers in the cluster may have very different arrival directions. In our opinion, the fact that the majority of EAS in the observed clusters are less “powerful” than an average shower witnesses in favour of this model since we do not need too strong magnetic fields. The duration of a cluster depends on the thickness of the LAP and the time during which the magnetic lens “works.”

Now let us ask ourselves how does it happen that the time series that represent the majority of EAS clusters (as well as the other data) are stochastic but some clusters seem to demonstrate signs of chaotic dynamics. At the moment, we cannot give a definite answer but can only conjecture that different factors may be involved in this phenomenon. Among them, one can mention a possibly fractal nature of the interstellar medium (see, e.g.,  and references therein), nonlinear mechanisms of particles acceleration, etc. A very interesting possibility for EAS time series to become chaotic suggests the fact that the large-scale fluctuations in the interplanetary magnetic field strength sometimes have fractal or multifractal structure . This phenomenon has been observed at different distances from Earth and for very different time scales. It is likely that in case that the “magnetic lens” discussed above has certain fractal or multifractal properties then it may lead to the appearance of chaotic dynamics in EAS time series.

Certainly, the presented model is only conceptual and oversimplified. Still we hope that it reflects the nature of processes that may lead to the appearance of EAS clusters, both “ordinary” and “possibly chaotic” ones.

Finally, it is interesting to compare our results with the conclusions of similar investigations performed by other research groups. In a considerable number of articles devoted to the nonlinear time series analysis, one can find a comprehensive investigation of EAS arrival times registered with the EAS-TOP array .

Basing on a detailed study of the available experimental data set and the results obtained with the underground muon monitor  the authors of this work made a conclusion that cosmic ray signals are all colour random noise, independently of the nature of the secondary particle and of the primary parent particle, but an existence of deterministic chaotic effects in cosmic ray time series cannot be completely excluded. It was also demonstrated in one of the following articles that an impact of background noise brings additional difficulties to the problem of distinguishing between chaotic and stochastic dynamics . In our opinion, these conclusions as a whole do not contradict our results, especially in view of the fact that signs of chaotic behaviour have only been observed for about 0.1% of EAS in the whole data set.

Besides this, a whole series of investigations devoted to the nonlinear analysis of EAS time series are carried out in Japan beginning from early 1990s at the experimental arrays that now constitute the LAAS network, see, e.g.,  and references therein. The authors of these investigations presented several dozens of events that demonstrate chaotic dynamics. More than this, it was conjectured that the observed dynamics may be due not only to the chaotic structure of the medium through which particles have traversed but also to the nature of the primary particles . Later on, there was suggested a model according to which chaotic events may be generated by cosmic rays that have a structure of a fractal wave arriving from a nonlinear accelerator like a supernova remnant . This model needs to be studied in details, but seems to be promising. Thus that the results obtained during our analysis do not contradict the conclusions of similar investigations performed at other EAS arrays.

The results presented above demonstrate that one can observe an unusual dynamics of EAS arrival times in the vicinity of certain clusters of EAS with the electron number of the order of $`10^5`$. While the overwhelming majority of samples in our data set are unambiguously stochastic, a number of samples in the vicinity of four EAS clusters demonstrate signs of chaotic dynamics. Still it is rather difficult to make a final conclusion on the nature of this phenomenon: Does it represent deterministic chaos or a special type of a stochastic process? In our opinion, the majority of the tests performed witness in favour of the first of these two alternatives. Nevertheless we must mention that our analysis may somehow suffer of the fact that the phenomena discussed above are only observed at comparatively short time scales with short samples while time series analysis usually prefers longer samples. In this connection, we recall that our investigation of EAS clusters has revealed an existence of “super clusters,” i.e., clusters that have duration more than 30 min and consist of hundreds of EAS. Our future plans include an analysis of these events.

It seems to be necessary to continue the work in this area and to involve some other methods of nonlinear time series analysis. There are a number of other nonlinear tools that may help to make a more definite conclusion about the nature of the observed phenomenon. Among them, one can recall space-time–separation and recurrence plots and the Lyapunov exponents. There are also a number of other measures of nonlinearity besides the one used above . Last but not least, special signal filtering techniques may be used in future to reduce effects of background noise on the dynamics of PC samples.

Finally, as we have already mentioned above, perhaps the biggest puzzle in connection with signs of chaos in EAS time series is their astrophysical nature. It is likely that clusters that produce signs of chaotic dynamics in the corresponding time series are similar to the upper part of an iceberg in a sense that they do not present the complete process but only the most pronounced part of it. We point out that for all PC samples discussed above, the value of the correlation dimension is comparatively small ($`D_2<3`$). Since this value gives a lower estimate for the number of degrees of freedom in the underlying process, it is possible that the structure of this process is not too complicated. Still it seems to be a great challenge to work out a good model that could explain chaotic dynamics in EAS arrival times.

We gratefully acknowledge numerous useful discussions with A. V. Igoshin, A. V. Shirokov, and V. P. Sulakov who have helped us a lot with the experimental data set. We also thank Thomas Schreiber for a very helpful communication and anonymous referees for stimulating comments. This work was done with financial support of the Federal Scientific-Technical Program “Research and design in the most important directions of science and techniques” for 2002–2006, contract No. 40.014.1.1.1110, and by Russian Foundation for Basic Research grant No. 02-02-16081. Only free, open source software was used for this investigation.


format: text

text:
# Magnetic enhancement of Co0.2Zn0.8Fe2O4 spinel oxide by mechanical milling

## I Introduction

In recent years, several research groups  are involved in the investigations of nanoparticle spinel oxides because of their potential applications in magnetic devices, in micro wave technology , in high density magnetic recording media , in magnetic fluids as drug carrier etc. . Various types of nanoparticle materials such as, metal: Fe, Co, Ni , metallic alloys: Fe-Cu , and metallic oxides: MnFe2O4  and ZnFe2O4 , are under current research activity. While metal and inter metallic nanoparticles suffer from stability problems in atmospheric condition, metallic oxides are highly stable under ambient conditions . Various factors such as, particle size distribution , inter-particle interactions , grain (core) and grain boundary (shell) structure  and metastable structure of the system  control the properties of nanoparticles. Some of the specific properties of the nanoparticles which are of interest are Quantum magnetic tunneling , various magnetic order like ferrimagnet/ferromagnet, spin glass/superparamagnet and spin canting effects .
The interesting aspect of magnetism in spinel oxides is that the magnetic order is strongly dependent on the competition between various superexchange type interactions i.e., JAB (A-O-B) and JBB (B-O-B), where A: tetrahedral (A) sites moments and B: octahedral (B) sites moments, O: is O2- ions . Certain amount of site disorder i.e., the cations redistribution between A and B sites is sufficient to change the super-exchange interactions in nanoparticles spinel . Hence, the magnetic properties of the nanoparticle spinels can be drastically different from their bulk counterpart. In this light, we investigated Co0.2Zn0.8Fe2O4 spinel oxide system.
The cation distribution of the bulk sample is (Zn$`{}_{0.8}{}^{}{}_{}{}^{2+}`$Fe$`{}_{0.2}{}^{}{}_{}{}^{+3}`$)A\[Co$`{}_{0.2}{}^{}{}_{}{}^{+2}`$Fe$`{}_{1.8}{}^{}{}_{}{}^{3+}`$\]BO4 ( A= tetrahedral sites, B= octahedral sites) . The A site is strongly diluted with respect to B site and B site moments contribute to different size clusters. The spins inside the clusters form canted structure due to short range ferrimagnetic interactions. The total magnetization M (= MBcos$`\theta `$-MA) of such type of canted system depends not only on the sublattice (site) magnetization MB and MA but also on the canting angle $`\theta `$ between the B site spins. The B site spin canting is strongly related with the amount of A site moments. If more number of A site moments exist in as milled Co0.2Zn0.8Fe2O4 sample, the canting angle between B site spins is expected to be reduced and a drastic change in cluster glass properties of bulk sample  is expected. Unfortunately, most of the reports in literature concern only on the long range order system, such as ZnFe2O4  and NiFe2O4 , and the magnitudes of A and B site moments but not on the short range interacting cluster glass system where change of spin canting inside the clusters plays an important role . In this presentation, we show that for as milled Co0.2Zn0.8Fe2O4 sample, the magnetization enhancement is related with the reduction of B site spin canting and an enhancement of ordering temperature occurs which is related with the cation redistribution and strain induced anisotropy.

## II Experimental

### A Sample preparation

Stoichiometric amounts of ZnO (99.998% from Johnson Matthey), Fe2O3 (99.998% from Johnson Matthey) and Co3O4 (99.5% from Fluka) oxides for the composition Co0.2Zn0.8Fe2O4 were mixed and ground for 2 hours. Then, the mixture was pelletized and sintered at 9500C for 12 hours and slowly cooled to room temperature. The pellet was again ground and pelletized and finally sintered at 10000C for 12 hours. The heat treatment was carried out under atmospheric condition with heating and cooling rate of 2-3 0C/minute. The XRD spectra of the as prepared bulk sample (S0) confirm formation of well crystalline cubic spinel phase. The as prepared sample was milled for 24 hours in a SPEX 8000 mixer/mill using a set of six balls, two of diameter 1/2 inch and four of diameter 1/4 inch with ball to powder weight ratio 5:1. The as milled sample (S1) was heat treated in air at 3000C for various durations, ranging from 2 hours to 575 hours, to study the temporal phase evolution. The 6 hours heat treated as milled samples are denoted as SN3 (for 3000C), SN6 (for 6000C) and SN10 (for 10000C). The heating and cooling rates were maintained at 2-30C/minute during the heat treatment process.

### B X-ray diffraction

The X ray diffraction (XRD) data (Fig. 1a) were taken using Philips PW1710 diffractometer with CuKα radiation in the 2$`\theta `$ range 100-900 with a step size 0.020. The XRD data of as prepared bulk sample show narrow crystalline lines with estimated particle size of a few micrometers. After 24 hours milling of the bulk sample, the XRD peak lines become significantly broad and no other phases are seen, except the cubic spinel structure. The broadening in XRD lines suggest nano structure or non-uniform micro strain introduced in the lattices during milling process . The width of $`<311>`$ XRD peak line gives average particle size $``$ 20 nm for the as milled sample (S1). When the as milled sample is heat treated, the XRD lines (Fig.1b-c) become sharp, which indicates recrystallization process in the system. The average particle size increases with annealing temperature as 23 nm, 28 nm and 62 nm for the samples SN3, SN6 and SN10, respectively. Eventhough the XRD spectrum of the SN10 sample is very similar to that of S0 sample, the average particle size $``$ 62 nm suggest that complete restoration of crystallinity close to S0 sample is yet to be reached.

### C X-ray Fluorescence spectra

We have determined the chemical composition of S1 sample by X-ray Fluorescence (XRF) technique. The Fluorescence spectrum was detected using ORTEC HPGe detector and finally recorded in a PC based multichhanel analyzer. The XRF spectrum (Fig. 2) shows the characteristic Kα and Kβ lines of Fe, Co and Zn elements. No extra lines have been found that correspond contamination from hardened steel balls with composition Fe74Cr18Ni8. The peak integrals of the Kα lines of all the elements are obtained using the standard peak fitting programme to get the relative percentage. This program uses the principle of Fundamental Parameter Method . The obtained composition is : Fe = 62.45%, Co = 7.17% and Zn = 30.38% which is very close to the expected values i.e. Fe = 63.53%, Co = 6.71% and Zn = 29.76% for Co0.2Zn0.8Fe2O4 system. The maximum error of the fitted values are within 5%.

### D Measurements

The low field dc magnetization (T = 20K to 320K, H= 10 Oe to 100 Oe) and ac susceptibility (T= 60K to 330K, hrms $``$ 1 Oe and frequency (f) = 337 Hz and 7 kHz) measurements have been performed using home made magnetometer . In ZFC condition, the sample has been cooled from room temperature to 20K in the absence of dc magnetic field, then the field was applied at 20K and magnetization data were recorded while increasing the temperature. In FC condition, the sample was cooled from room temperature to 20K in presence of dc magnetic field (same magnitude which was applied during ZFC measurement) and the magnetization data were recorded with increasing temperature and keeping the field on. High field magnetization, hysteresis experiments were performed using VSM magnetometer in the fields upto $`\pm `$12 Tesla in the temperature range 10K to 300K.
Mössbauer spectra were recorded at 300K without applying external magnetic field, using a constant acceleration spectrometer in transmission geometry mode. The spectra were recorded using a 5 mCi 57Co in Rh matrix source. The hyperfine magnetic field (HMF) distribution at the 57Fe nuclei was evaluated from the Mössbauer spectra using the method of Le Caër and Dubois . In this model a linear relationship between HMF (H) and isomer shift (IS) is assumed in the form IS = aH + b, where a and b are the fitting parameters to get a minimum $`\chi ^2`$. The isomer shift (IS) wes calculated with respect to 57Fe metal.

## III Results and Discussion

### A AC susceptibility

Fig. 3 shows the real ($`\chi ^{}`$) and imaginary ($`\chi ^{\prime \prime }`$) components of ac susceptibility for sample S1. It shows a broad $`\chi ^{}`$ maximum at TB $``$ 320K (for hrms $``$ 1 Oe, f= 337 Hz) and the corresponding $`\chi ^{\prime \prime }`$ maximum at $``$ 290K. The broadness of the ac susceptibility maxima suggest a cluster size distribution in the as milled sample S1. The temperature where $`\chi ^{}`$ shows maximum is defined as the average blocking temperature (TB) of the clusters. From the frequency dependence of ac susceptibility measurements (Fig. 3a), a positive frequency shift of TB is observed in the $`\chi ^{}`$ as well as in the $`\chi ^{\prime \prime }`$ maximum in the measurement frequency range. This is a characteristic of superparamagnetic blocking of the ferromagnetic clusters in different metastable states . In the low temperature region, the $`\chi ^{}`$ data show a small shoulder at $``$ 100K (designated as TK) and $`\chi ^{\prime \prime }`$ data show (Fig. 3b) increasing tendency below TK. This type of behaviour of TK can be attributed due to the disordered surface or grain boundary spins in a nanoparticle system .
In order to check the effect of heat treatment on magnetic properties, we have carried out the ac susceptibility measurements for the SN10 sample. This heat treated sample shows mixed magnetic phases. The $`\chi ^{}`$ and $`\chi ^{\prime \prime }`$ data of SN10 sample (Fig. 4) show two magnetic ordering at Tm1 $``$ 70K and Tm2 $``$ 280K respectively. The frequency shift of Tf follows Vogel-Fulcher law

$$f=f_0exp^{E_a/(T_{m1}T_0)}$$
(1)
with characteristic frequency $`f_0`$ $``$ 1010 Hz, activation energy Ea $``$ 379 eV and constant T0 $``$ 47K. The corresponding frequency shift per decade of f ($`\mathrm{\Delta }T_f`$/T$`{}_{f}{}^{}\mathrm{\Delta }log(f)`$) is $``$ 0.05, which characterizes spin glass like ordering at Tm1 . The ordering temperature at Tm2 does not show significant frequency shift (not shown in Fig. 4), suggest ferrimagnetic ordering temperature at Tm2 $``$ 280K. The spin glass transition temperature Tm1 $``$ 70K and ferrimagnetic ordering temperature at Tm2 $``$ 280K of SN10 sample are slightly different with respect to the cluster spin glass freezing temperature Tm1 $``$ 110K and short range ferrimagnetic ordering temperature Tm $``$ 260K for bulk sample . But it is clear that the mixture of two magnetic phases (spin glass/cluster spin glass plus ferrimagnetic phase) of SN10 sample is very similar to the bulk sample S0.

### B DC magnetization

The ZFC magnetization of S1 sample (Fig. 5) shows a broad maximum at the cluster blocking temperature TB $``$ 320K (H $``$ 30 Oe) with a thermomagnetic irreversibility between ZFC and FC magnetization below TB. The decrease of ZFC magnetization below TB is due to the blocking of ferromagnetic clusters in different metastable states, whereas the FC magnetization increases below TB due to the orientation of these ferromagnetic clusters in metastable states which give rise to more magnetic contribution . It is observed (Fig.5a and Fig. 5b) that the TB is highly applied field dependent (i.e., TB $``$ 320K, 280K and 50K for applied field 30 Oe, 100 Oe and 1 Tesla, respectively) as expected for superparamagnetic blocking of the clusters . The superparamagnetic behaviour of the as milled sample (S1) can be further identified from the temperature dependence of field cooled thermoremanent magnetization (TRM) data (Fig. 5b inset) that the TRM value reduces to zero at T$``$ TB. The zero TRM value indicate that the effective inter-cluster interactions are also negligible and the clusters behave as a non-interacting small particle above TB $``$ 320K.
An interesting feature is observed when we compare the magnetization data at 100 Oe of as prepared bulk sample, as milled sample and heat treated samples (Fig. 6). The magnetization is enhanced when the bulk sample is milled for 24 hours. On heat treating the as milled sample, it is observed that magnetization further increases for SN3 sample. Then magnetization decreases for the SN6 sample. If the as milled sample is heat treated at higher temperature (SN10), the magnetization increases. But compared to the bulk sample S0, SN10 sample has a higher magnetization value at low temperature and lower value at temperatures $`>`$ 200K. Similar non-equilibrium magnetic behaviour as a function of annealing temperature has been observed for NiFe2O4 . The superparamagnetic blocking temperature TB (indicated in Fig. 6) decreases to 270K and 120K for SN3 and SN6, respectively, in comparison with TB $``$ 280K for S1 sample at H= 100 Oe. The 100 Oe magnetization data of the SN10 sample indicate a mixture of two magnetic ordering at Tm1 $``$ 70K and Tm2 $``$ 225K, respectively, very similar to the cluster spin freezing temperature at Tm1 $``$ 100K and short range ferrimagnetic ordering temperature at Tm2 $``$ 230K, respectively, for bulk sample .

### C Hysteresis

The field dependence of magnetization data (Fig. 7) under ZFC condition is similar to that of ferromagnetic isotherms but the magnetization lacks saturation even upto 12 Tesla at any temperature. This, we attribute to the spin canting effects at grain boundary  or superparamagnetic contribution of nanoparticles . We note from Fig. 7a (inset) that the S1 sample has a better ferromagnetic behaviour in terms of saturation with respect to magnetic field which indicates that the strong spin canting behaviour of S0 sample  has decreased for the nanoparticle S1 sample. Importantly, we also note (Fig.7a, inset) that the saturation magnetization of nanoparticle sample S1 is significantly increased at room temperature with respect to the bulk sample S0. However, the higher value of M for S0 at 10K for H $`>`$ 4T compared to S1 can be understood by assuming two types of magnetic interactions in the system. One is ferromagnetic and second one is antiferromagnetic. In S0 sample the B site antiferromagnetic interactions are dominant which causes spin canting in B site and shows non-saturation in magnetization . The reduction of antiferromagnetic interactions between the B site moments gives rise the ’better’ ferromagnetic behaviour, yet having a reduced saturation moment at 10K, for S1 sample. It is observed from the hysteresis loop (Fig. 7b) that isothermal remanent magnetization (MR), coercive field (HC) where MR is zero, and irreversible field (Hirr), where the hysteresis loop closes, decrease with temperature. There is no hysteresis loop at T $``$ 290K. This suggest that the thermal energy close to the blocking temperature ($`T_B`$) is sufficient to reduce the coercive field to zero value . The saturation magnetization (MS) values, obtained from M vs 1/H plot, are small in comparison with the bulk sample as seen in nanocrystalline NiFe2O4 . This is possible if the B site intra-cluster spin canting decreases or if there is a redistribution of cations between A and B sites when the bulk sample is mechanically milled .

### D Mössbauer Spectroscopy

Fig. 8 shows the room temperature Mössbauer spectra in absence of any external magnetic field for the samples (S0), S1, SN3, SN6 and SN10. The spectrum of S0 sample consists of a Lorentzian doublet arising from the Fe3+ ions at octahedral (B) site and a Lorentzian singlet arising from the Fe3+ ions at tetrahedral (A) site. The most probable cations distribution of bulk sample (S0) obtained is (Zn0.8Fe0.2)A\[Co0.2Fe1.8\]BO4. The isomer shift (IS) and quadrupole splitting (QS) of Fe3+ ions at \[B\] site are $``$ +0.29 mm/sec and 0.35 mm/sec respectively. The IS value of Fe3+ ions at (A) site is $``$ -0.056 mm/sec . The Mössbauer spectrum of the bulk sample therefore confirm that Zn2+ ions occupy A site, Co2+ ions prefer B site and Fe3+ ions prefers both A and B sites. The Mössbauer spectrum of the S1 sample clearly indicates a hyperfine magnetic field splitting in addition to a central paramagnetic doublet. This spectrum of S1 sample represents the appearance of spontaneous magnetization of ferromagnetic clusters mixed with superparamagnetic fluctuation effect due to nano meter size of the particles . As the annealing temperature increases, the hyperfine magnetic field splitting decreases for SN3 sample and the SN6 and SN10 samples show paramagnetic spectra. The Mössbauer parameters of SN10 sample (IS= +0.27 mm/sec and QS = +0.43 mm/sec for B site Fe3+ ions and the IS = + 0.25 mm/sec of A site Fe3+ ions) differ from that of the as prepared bulk sample which suggests that duration of heat treatment ($``$ 6 hours) is not enough to achieve 100% equilibrium state of the bulk samples. These results confirm that the enhancement of magnetic ordering and magnetization observed in S1 sample is intrinsic property of the material. These also further indirectly confirm the hypothesis of the presence of non-equilibrium cation distribution in the as milled sample.
To check the recovery of cation distribution from non-equilibrium state to equilibrium state, we recorded the Mössbauer spectrum of S1 sample as function of duration of heat treatment at 3000C. Fig. 9 shows the Mössbauer spectra and corresponding hyperfine magnetic field (HMF) distribution (p(H)). The presence of multiple hyperfine fields at A and B site Fe3+ ions can be expected due to various super transferred hyperfine fields (STH) from neighbouring ions. The HMF acting on the B site Fe3+ ions is due to non-uniform environment, consisting of different number of Zn2+ and Fe3+ ions at the nearest neighbour A sites. The HMF experienced by the A site Fe3+ ions is due to nearest neighbour B site configuration . It is clear from the p(H) distribution that the changes brought about by 2 hours heat treatment is not very significant but the peak intensities of hyperfine field components are decreasing and intensities of paramagnetic components are increasing for annealing time $``$ 16 hours. Indirectly, it gives the information that more and more number of Zn2+ ions stabilizing in A site with the compromise of more Fe3+ ions in B site. This reduces the reduction of super transferred hyperfine magnetic field via inter-sublattice superexchange interactions . Fig. 10 shows the temporal evolution of the average HMF for 3000C heat treated sample. The solid line represent the fit data using a function of time (t) as F(t) = a \+ b(1-exp(t/t0)) where a, b and t0 are constants. F(t) represents average HMF as a function of time. The fit gives a measure of the reaction kinetics through atomic diffusion as the sample gets annealed.

## IV Summary and Conclusions

Usually nano materials lead to the decrease in magnetization and ordering temperature as the particle size decreases . This is already observed in Co0.2Zn0.8Fe2O4 with particle size $``$ 6nm to 70 nm, prepared by coprecipitation method . However, the bulk Co0.2Zn0.8Fe2O4 sample with short range ferrimagnetic ordering temperature $``$ 260K and paramagnetic at 300K  after 24 hours mechanical milling gives nanoparticles $``$ 20 nm which shows ferromagnetic cluster blocking state below TB $``$ 320K with an enhancement of magnetization. Because of the weak inter-sublattice interactions JAB (due to low A site Fe3+ population), the B site spins form canted structure in bulk Co0.2Zn0.8Fe2O4 sample. The stable cation distribution of the bulk sample is destroyed due to mechanical milling and a non-equilibrium disordered state, i.e. the increase of A site Fe3+ population enhances the inter-sublattice exchange interactions Fe$`{}_{A}{}^{}{}_{}{}^{3+}`$-O2--Fe$`{}_{B}{}^{}{}_{}{}^{3+}`$ (JAB) and change the B sublattice superexchange interactions Fe$`{}_{B}{}^{}{}_{}{}^{3+}`$-O2--Fe$`{}_{B}{}^{}{}_{}{}^{3+}`$. This makes the clusters more ferromagnetic with lower spin canting. The reduction of spin canting inside the clusters enhances the magnetization. The superparamagnetic fluctuation effect is expected due to the nano size clusters (particles). The cations redistribution surely is the main factor which controls the magnetic properties of nanoparticle system but distortion of spherical shape and strain induced anisotropy energy due to mechanical milling can give rise to a preferential orientation of the spins inside the clusters.
The magnetic reversibility of the heat treated (SN10) milled sample towards the bulk sample (may not 100%, due to different micro structural parameters like particle size, cluster size distribution, spin canting angle inside the clusters) conclusively show that enhancement of magnetization and ordering temperature of the as milled samples are intrinsic property of Co0.2Zn0.8Fe2O4 spinel oxide.
Acknowledgement: One of the authors RNB thanks Council of Scientific and Industrial Research (CSIR, India) for providing fellowship \[F.No.9/489(30)/98-EMR-I \].

Figure Caption

Fig.1 XRD spectra for Co0.2Zn0.8Fe2O4 spinel oxide heated at different temperatures. For the notation S0, S1, SN3, SN6 and SN10 see in the text.
Fig.2 X-ray Fluorescence Spectrum for the as milled sample. the arrow indicates the corresponding metals characteristics Kα and Kβ lines.
Fig.3 Ac susceptibity $`\chi ^{}`$ (in a) and $`\chi ^{\prime \prime }`$ (in b) data for as milled sample (S1). TB and TK represent blocking and spin canting temperature, respectively.
Fig.4 Ac susceptibity $`\chi ^{}`$ (in a) and $`\chi ^{\prime \prime }`$ (in b) for the 10000C-6 hour heat treated sample (SN10). Fig.a inset shows ac susceptibility data at 1 Oe ac field, f = 337 Hz and ln f vs 1/Tm1 plot for the same sample.
Fig.5 Dc magnetization vs temperature for 24 hrs milled Co0.2Zn0.8Fe2O4 sample. Fig. b inset shows the temperature dependence of field cooled (FC) (cooled in a field of 60 Oe) thermoremanent magnetization (TRM).
Fig.6 Temperature dependence of dc magnetization data for different heat treated samples at 100 Oe field.
Fig.7 a) M vs H plot at different temperatures for S1 sample and inset Fig. compare the 10K and 290K data of S0 and S1 samples. b) Hysteresis loopes at different temperatures for S1 sample. For MR, HC and Hirr see text.
Fig.8 Room temperature Mössbauer spectra recorded in absence of field for as prepared bulk (S0) and heat treated as milled 1, SN3, SN6 and SN10 sample.
Fig.9 Room temperature Mössbauer spectra for heat treated Co0.2Zn0.8Fe2O4 sample at 3000C for different time (left) and corresponding hyperfine field distribution (right).
Fig.10 Annealing time (at 3000C) dependence of average hyperfine field value for as milled Co0.2Zn0.8Fe2O4 sample.


text:
# The Sunyaev-Zel’dovich effect revisited

## 1 Introduction

As it is well known today, the cosmic microwave background radiation (CMBR) that fills the universe exhibits an almost perfect blackbody Planckian spectrum corresponding to a temperature of 2.726 K. A little over thirty years ago, it was predicted by the astrophysicists R.A Sunyaev and Ya. B. Zel’dovich that this spectrum could be distorted when photons from this radiation would penetrate large structures such as the hot intracluster gas now known to exist in the universe . This distortion would arise from the interaction between the photons and electrons which constitute the plasma responsible for an important main contribution to the total mass contained in the cluster. Such distortion is only a very small effect changing the brightness of the spectrum by a figure of the order of 0.1 percent. This effect is now called the Sunyaev-Zel’dovich effect and its detection is at present a relatively feasible task due to the modern observational techniques available. Its main interest lies on the fact that it provides information to determine important cosmological parameters such as Hubble’s constant and the baryonic density. Since all these and other important facts are readily available in the somewhat broad literature on the subject -, we will not pursue any further details of its cosmological implications.

In this work we are mainly concerned with the nature of the various interpretations that have been provided to explain the effect and to propose another one which we believe, is much simpler to grasp than the others, and offers a much more direct way to account for the relativistic corrections, as well as the influence of $`z`$, the redshift factor, in the equations for the distorted spectrum.

The main idea behind the method we want to discuss is that photons incident on the hot electron gas of the plasma change their frequency by an absortion-emission process, so that the intensity of the spectral line is changed through the Doppler effect.

The overall effect is due to those photons that happen to be captured by an electron, an event which is in general unlikely to occur specially in an optically thin gas. Thus, an electron moving with a given thermal velocity emits (scatters) a photon with a certain incoming frequency $`\nu _o`$ and outgoing frequency $`\nu `$. The line breath of this process is readily calculated from elementary kinetic theory taking into account that the media in which the process takes place has an optical depth directly related to the usual Compton parameter $`y`$. When the resulting expression is convoluted with the incoming flux of photons obtained from Planck’s distribution, one easily obtains the disturbed spectra. The results derived from the final equation are shown to be in agreement with those obtained from the observational data. To present these results we have divided the paper as follows. In section II a brief review of the previous methods developed to explain the SZ effect is given. This will allow a thorough comparison with our proposal. In Section III we develop our ideas and give the comparison between the mathematical results with both, those derived before and the observational ones. Section IV is left for some concluding remarks and future directions of this work.

## 2 The Diffusion-Scattering Picture of the SZ effect.

As it was clearly emphasized by the authors of this discovery in their early publications   as well as by other authors, the distortion in the CMBR spectrum by the interaction of the photons with the electrons in the hot plasma filling the intergalactic space is due to the diffusion of the photons in the plasma which, when colliding with the isotropic distribution of a non-relativistic electron gas, generates a random walk. The kinetic equation used to describe this process was one first derived by Kompaneets back in 1956 -. For the particular case of interest here, when the electron temperature $`Te10^8K`$ is much larger than that of the radiation, $`TT_{Rad}(2.726K)`$, the kinetic equation reads

$$\frac{dN}{dy}=\nu ^2\frac{d^2N}{d\nu ^2}+4\nu \frac{dN}{d\nu }$$
(1)
Here $`N`$ is the Bose factor, $`N=(e^x1)^1`$ , $`x=\frac{h\nu }{kT}`$, $`\nu `$ is the frequency, $`h`$ is Planck’s constant, $`k`$ Boltzmann’s constant and $`y`$ the “Compton parameter” given by,

$$y=\frac{kTe}{m_ec^2}\sigma _Tn_ec𝑑t=\frac{kTe}{m_ec^2}\tau $$
(2)
$`m_e`$ being the electron mass, $`c`$ the velocity of light, $`\sigma _T`$ the Thomson’s scattering cross section and $`n_e`$ the electron number density. The integral in Eq. (2) is usually referred to as the “optical depth, $`\tau `$”, measuring essentially how far in the plasma can a photon travel before being captured (scattered) by an electron. Without stressing the important consequences of Eq. (1) readily available in many books and articles -, we only want to state here, for the sake of future comparison, that since the observed value of $`N_o(\nu )`$ is almost the same as its equilibrium value $`N_{eq}(\nu )`$, to a first approximation one can easily show that

$$\frac{\delta N}{N}\frac{N_o(\nu )N_{eq}(\nu )}{N_{eq}(\nu )}=y\left[\frac{x^2e^x(e^x+1)}{(e^x1)^2}\frac{4xe^x}{e^x1}\right]$$
(3)
implying that, in the Rayleigh-Jeans region $`\left(x<<1\right)`$,

$$\frac{\delta N}{N}2y$$
(4)
and in the Wien limit $`\left(x>>1\right)`$ :

$$\frac{\delta N}{N}x^2y$$
(5)
If we now call $`I_o(\nu )`$ the corresponding radiation flux for frequency $`\nu `$, defined as

$$I_o(\nu )=\frac{c}{4\pi }U_v(T)$$
(6)
where $`U_v(T)=\frac{8\pi h\nu ^3}{c^3}N_{eq}(\nu )`$ is the energy density for frequency $`\nu `$ and temperature $`T`$ , and noticing that

$$\frac{\delta T}{T}=\left(\frac{(\mathrm{ln}I(\nu ))}{(\mathrm{ln}T)}\right)\left(\frac{\delta I}{I}\right)$$
(7)
we have that the change in the background brightness temperature is given, in the two limits, by

$$\frac{\delta T}{T}2y,x<<1$$
(8)
and

$$\frac{\delta T}{T}xy,x>>1$$
(9)
showing a decrease in the low frequency limit and an increase in the high frequency one. Finally, we remind the reader that the curves extracted from Eq. (3) for reasonable values of the parameter $`\mathrm{`}\mathrm{`}y\mathrm{"}`$ show a very good agreement with the observational data.

Nevertheless, many authors, including Sunyaev and Zel’dovich themselves, were very reluctant in accepting a diffusive mechanism as the underlying phenomena responsible for the spectrum distortion. The difficulties of using diffusion mechanisms to study the migration of photons in turbid media, specially thin media, have been thoroughly underlined in the literature --. Several alternatives were discussed in a review article in 1980  and a model was set forth by Sunyaev in the same year  based on the idea that Compton scattering between photons and electrons induce a change in their frequency through the Doppler effect. Why this line of thought has not been pursued, or at least, not widely recognized, is hard to understand. Two years ago, one of us (ASV)  reconsidered the single scattering approach to study the SZ effect. The central idea in that paper is that in a dilute gas, the scattering law is given by what in statistical physics is known as the  *dynamic structure factor*, denoted by $`S(k,\nu )`$ where $`k=\frac{2\pi }{\lambda }`$, and $`\lambda `$ is the wavelength. In such a system, this turns out to be proportional to $`\mathrm{exp}(\frac{\nu ^2}{w^2})`$ where $`w`$ is the broadening of the spectral line given by

$$w=\frac{2}{c}(\frac{2kTe}{m})^{1/2}\nu $$
(10)
If one computes the distorted spectrum through the convolution integral

$$I\left(\nu \right)=_0^{\mathrm{}}I_o\left(\overline{\nu }\right)S(k,\overline{\nu }(1ay)\nu )𝑑\overline{\nu }$$
(11)
where the corresponding frequency shift $`\frac{\delta \nu }{\nu }`$ has been introduced through Eqs. (8-9) $`(a=2+x)`$, one gets a good agreement with the observational data. This is exhibited in Fig. (1) for $`y=10^5`$. This result is interesting from, at least, two facts. One, that such a simple procedure is in agreement with the diffusive picture. This poses interesting mathematical questions which will be analyzed elsewhere, specially since the exact solution to Eq. (1) is known (see Eq. (A-8) ref. ). The other one arises from the fact that this is what triggered the idea of reanalyzing the SZ effect using elementary arguments of statistical mechanics and constitutes the core of this paper to be presented next.

## 3 The Doppler Effect Approach.

We begin this section by simply reminding the reader that if an atom in an ideal gas moving say with speed $`u_x`$ in the $`x`$ direction emits light of frequency $`\nu _o`$ at some initial speed $`u_x(0)`$, the intensity of the spectral line $`I\left(\nu \right)`$ is given by

$$\frac{I\left(\nu \right)}{I_o\left(\nu \right)}=\mathrm{exp}\left[\frac{mc^2}{2kT}\left(\frac{\nu _o\nu }{\nu }\right)^2\right]$$
(12)
Eq. (12) follows directly from the fact the velocity distribution function an ideal gas is Maxwellian and that $`u_x`$, $`u_x(0)`$ and $`\nu `$ are related through the Doppler effect. For the case of a beam of photons of intensity $`I_o\left(\nu \right)`$ incident on a hot electron gas regarded as an ideal gas in equilibrium at a temperature $`T_e`$ the full distorted spectrum may be computed from the convolution integral given by

$$I\left(\nu \right)=\frac{1}{\sqrt{\pi }W\left(\nu \right)}_0^{\mathrm{}}I_o\left(\overline{\nu }\right)\mathrm{exp}\left[\left(\frac{\overline{\nu }f(y)\nu }{W\left(\nu \right)}\right)^2\right]𝑑\overline{\nu }$$
(13)
which defines the joint probability of finding an electron scattering a photon with incoming frequency $`\overline{\nu }`$ and outgoing frequency $`\nu `$ multiplied by the total number of incoming photons with frequency $`\overline{\nu }`$. $`\frac{1}{\sqrt{\pi }}W\left(\nu \right)`$ is the normalizing factor of the Gaussian function for $`f(y)=1`$, $`W\left(\nu \right)`$ is the width of the spectral line at frequency $`\nu `$ and its squared value follows from Eq.(10)

$$W^2\left(\nu \right)=\frac{4kT_e}{m_ec^2}\tau \nu ^2=4yv^2$$
(14)
where $`\tau `$ is the optical depth whose presence in Eq. (14) will be discussed later. The function $`f(y)`$ multiplying $`\nu `$ in Eq. (13) is given by $`f(y)=1+ay`$ where $`a=2`$ in the Rayleigh-Jeans limit and $`a=xy`$ in the Wien’s limit, according to Eqs. (8-9) and the fact that $`\frac{\mathrm{\Delta }\nu }{\nu }=\frac{\mathrm{\Delta }T}{T}`$ for photons. Eq.(13) is the central object of this paper so it deserves a rather detailed examination. In the first place it is worth noticing that $`I_o(\nu )`$, the incoming flux,is defined in Eq. (6). Secondly, it is important to examine the behavior of the full distorted spectrum in both the short and high frequency limits. In the low frequency limit, the Rayleigh-Jeans limit $`I_o(\nu )`$ $`=\frac{2kT\nu ^2}{c^2}`$, so that performing the integration with $`a=2`$ and noticing that $`y`$ is a very small number, one arrives at the result

$$\frac{\delta I}{I}\frac{I(\nu )I_o(\nu )}{I_o(\nu )}=2y,x<<1$$
(15)
(15) is in complete agreement with the value obtained using Eq.( 3), the photon diffusion equation. In the high frequency limit where $`a=xy`$ and $`I_o(\nu )`$ $`=\frac{2h\nu ^3}{c^2}e^x`$, a slightly more tedious sequence of integrations leads also to a result at grips with the diffusion equation, namely

$$\frac{\delta I}{I}=x^2y,x>>1$$
(16)
Why both asymptotic results, the ones obtained with the diffusion equation and those obtained from Eq. (13) agree so well, still puzzles us. At this moment we will simply think of them as a mathematical coincidence. Nevertheless, it should be stressed that in his 1980 paper, Sunyaev  reached rather similar conclusions although with a much more sophisticated method, and less numerical accuracy for the full distortion curves. The distorted spectrum for the CMBR radiation may be easily obtained by numerical integration of Eq. (13) once the optical depth is fixed, $`y`$ is determined through Eq.(2) and $`a=2+x`$. The intergalactic gas cloud in clusters of galaxies has an optical depth $`\tau 10^2`$ . To compare our results with the curves obtained from the traditional approach, we have plotted $`\delta I\left(\nu \right)`$ for several temperatures, typical of the hot intra-cluster gas, in the non-relativistic range, in Figs. (2) and (3).

From the results obtained one appreciates the rather encouraging agreement between the observational data and the theoretical results obtained with the three methods, the diffusion equation, the structure factor or scattering law approach, and the Doppler effect. This, in our opinion is rather rewarding and some efforts are in progress to prove the mathematical equivalence of the three approaches. From the physical point of view, and for reasons already given by many authors, we believe that the scattering Doppler effect picture does correspond more with reality, specially for reasons that will become clear in the last section.

## 4 Discussion of the results

Eq. (13), the main result of this work hardly needs a more detailed explanation, it is a direct result of elementary statistical mechanics, except of course for the factor $`\tau `$ which appears in the broadening of the spectral line. That this broadening must somehow depend on $`\tau `$ is clear. As argued by many authors in the case of clusters of galaxies, the plasma is optically thin, so most of the photons are not scattered, and how many of them are must depend on the optical depth. Now, why in particular Eq. (14) is valid remains to be rigorously shown. In our case it came as a mere accident since in ref.  the mass taken in $`w`$ to perform the calculations was that of a proton, which turns out to be of the same order of magnitude as $`\frac{m_e}{\tau }`$ for $`y=10^5`$. There is in fact one way of understanding this puzzle. In a pure diffusive (Fickian) process described by a Gaussian function one knows that the squared value of a single line width grows as $`t`$. Looking at the original equation of Kompaneets and thinking only in its diffusive terms, one sees that the variable time is replaced by the Compton parameter ”$`y`$”, so that one should expect that, to a first approximation, the width of the curve grows as $`\sqrt{y}`$. Hence, inside the gas, the effective frequency at which the photon propagates through the gas is, according to Eq. (14), proportional to $`\sqrt{y}\times \nu `$. This argument, which is based on ideas pertinent to the diffusion mechanism must be extended to the Doppler picture, but this has not yet been accomplished.

The two main advantages of Eq. (13) are, firstly, that its generalization to the relativistic case is straightforward and indeed, since in relativity theory one has both the transverse and parallel Doppler effects , one can easily study what their influence is, if any, in the full distorted spectrum. For the parallel effect the calculations carried out so far agrees reasonably well with those reported by Rephaeli and other authors  and will be the subject of a forthcoming publication. The case of the transverse effect is under study. Secondly, in all the work we have reported here, the cluster of galaxies, the source of the scattering processes of the CMBR and the electrons, has been considered a system at rest. One can now extend the calculations assuming that the cluster recedes to infinity according to Hubble’s law and introduce the redshift factor into the formalism. The resulting prediction could be compared with the observational data and extract valuable information on such an important cosmological parameter. Work in this direction is also in progress. Concluding, we have the strong conviction that the Doppler effect approach to study the thermal SZ effect is not only very clear from the physical point of view, it also contains ingredients which are promising for further studies more advantageous than those offered by the diffusive model.