format: text

text:
# Q-balls in Underground Experiments

## 1 INTRODUCTION

Dark Matter (DM) is one of the most intriguing problems in particle physics and cosmology. Several types of stable particles hypothesized in theories beyond the Standard Model of particle physics have been considered as candidates for DM. One example of such particles is the lightest supersymmetric particle (LSP) coming from a supersymmetric extension of the Standard Model .

In theories where scalar fields, carry a conserved global quantum number, $`Q`$, there may exist non-topological solitons which are stabilized by global charge conservation. These particles are spherically symmetric and for large values of $`Q`$ their masses and volumes grow linearly with $`Q`$. Thus they act like homogenous balls of ordinary matter, with $`Q`$ playing the role of particle number; Coleman called this type of matter Q-balls .

The conditions for the existence of absolutely stable Q-balls are satisfied in supersymmetric theories with low energy supersymmetry breaking. According to  abelian non-topological solitons with Baryon and/or Lepton quantum numbers naturally appear in the spectrum of the Minimal Supersymmetric Standard Model. The role of conserved quantum number is played by the baryon number. The same reasoning applies to sleptons for the lepton number and also to scalar Higgs particles. Q-balls can thus be considered like coherent states of squarks, sleptons and Higgs fields. Under certain assumptions about the internal self interaction of these particles and field the Q-balls are absolutely stable .

In this note we recall the main physical and astrophysical properties of these particles, the interaction of Q-balls with matter and their energy losses in matter, and the possibility of traversing the Earth to reach the MACRO detector . We neglect the possibility of :
$`(i)`$ electromagnetic radiation emitted by Q-balls of high $`\beta `$.
$`(ii)`$ strong interaction of Q-balls in the upper atmosphere capable of destroying the Q-ball ; this last point deserves further investigations.

## 2 PROPERTIES OF Q-BALLS

Q-balls could have been produced in the Early Universe, and could contribute to the DM. Several mechanisms could have lead to the formation of Q-balls in the Early Universe. They may have been created in the course of a phase transition, which is sometime called “solitogenesis”, or they could have been produced via fusion in processes reminiscent of the big bang nucleosynthesis, wich have been called “ solitosynthesis”; small Q-balls can be pair-produced in high energy collisions .

The astrophysical consequenses of Q-balls in many ways resemble those of strange quark matter, “nuclearites”; the peculiarity of Q-balls is that their mass grows as $`Q^{3/4}`$, while for nuclearites the mass grows linearly with baryon number .

The Q-ball mass and size are related to its baryon number . For a supersymetic potential $`U(\varphi )M_S^4=constant`$ for large scalar $`\varphi `$, the Q-ball mass $`M`$ and radius $`R`$ are given by .

$$M=\frac{4\pi \sqrt{2}}{3}M_SQ^{3/4}$$
(1)
$$R=\frac{1}{\sqrt{2}}M_S^1Q^{1/4}$$
(2)
The parameter $`M_S`$ is the energy scale of the SUSY breaking symmetry. A stability condition was found in ref. : the Q-ball mass $`M`$ is related to the nucleon mass $`M_N`$ by

$$MQM_N$$
(3)
From Eq. 1 and Eq. 3 one has the stability constraint :

$$Q1.6\times 10^{15}\left(\frac{M_S}{1TeV}\right)^4$$
(4)
In Fig. 1 the allowed region for stable Q-balls is indicated.

Q-balls are expected to concentrate in galactic halos and to move with the typical galactic velocity $`v=\beta c10^3c`$. Assuming that Q-balls constitute the cold galactic dark matter with $`\rho _{DM}0.3GeV/cm^3`$, their number density is

$$N_Q\frac{\rho _{DM}}{M}5\times 10^5Q^{3/4}\left(\frac{1TeV}{M_S}\right)cm^3$$
(5)
The corresponding flux is

$$\varphi 10^2Q^{3/4}\left(\frac{M_S}{1TeV}\right)^1cm^2s^1sr^1$$
(6)
If we assume that $`\beta _{Qball}10^3`$, $`\rho _{DM}0.3GeV/cm^3`$, $`M_S=1TeV`$, the Q-ball flux cannot be greater that $`4\times \frac{10^{19}}{M_Q}cm^2s^1sr^1`$ ($`M_Q`$ in g). Q-balls are expected to be distributed uniformly in our part of the galaxy and there should not be enhancements in the solar system, as for example a cloud of Q-balls around the sun.

Q-balls can be classified in two classes: SECS (Supersymmetric Electrically Charged Solitons) and SENS (Supersymmetric Electrically Neutral Solitons). The interaction of Q-balls with matter and their detection differ drastically for SECS or SENS.

## 3 INTERACTION WITH MATTER

### 3.1 Interaction with matter of Q-balls type SECS

SECS are Q-balls with a net positive electric charge in the interior. The charge of SECS originates from the unequal rate of absorption in the condensate of quarks (squarks) and electrons (selectrons). This positive electric charge is neutralized by a surrounding cloud of electrons. The positive charge interacts via elastic or quasi elastic collisions. The positive electric charge can vary from one to several tens; the cross section is the Bohr cross section for Q-ball hydrogen-interaction

$$\sigma =\pi a_0^210^{16}cm^2$$
(7)
Where $`a_0`$ is the Bohr radius. The formula is valid for $`Ra_0`$, which happens for $`Q2.7\times \left(\frac{M_S}{1TeV}\right)^4`$.

The main energy losses  of SECS passing throuth matter with velocities in the range $`10^4<\beta <10^2`$ is due to two contributions: the energy losses due to $`(i)`$ the interaction of the SECS core with the nuclei (nuclear contribution) and, $`(ii)`$ with the electrons of the traversed medium (electronic contribution). The total energy loss is the sum of the two contributions.

SECS could cause the catalysis of proton decay, but only if they are large and have large velocities . The possibility that SECS can cause the catalysis of proton decay does not concern our range of interest for velocities, masses and radii of SECS.

Electronic losses of SECS: The electronic contribution to the energy loss of SECS is calculated with the following formula

$$\frac{dE}{dx}=\frac{8\pi a_0e^2\beta }{\alpha }\frac{Z_1^{7/6}N_e}{(Z_1^{2/3}+Z_2^{2/3})^{3/2}}forZ_11$$
(8)
where $`\alpha `$ is the fine structure constant, $`\beta =v/c`$, $`Z_1`$ is the positive core charge of SECS, $`Z_2`$ is the atomic number of the medium and $`N_e`$ is the density of electrons in the medium. Electronic losses dominate for $`\beta >10^4`$.

Nuclear losses of SECS: The nuclear contribution to the energy loss of SECS is due to the interaction of the SECS positive core with the nuclei of the medium and it is given by

$$\frac{dE}{dx}=\frac{\pi a^2\gamma NE}{ϵ}S_n(ϵ)$$
(9)
where

$$S_n(ϵ)\frac{0.56Log(1.2ϵ)}{1.2ϵ(1.2ϵ)^{0.63}},ϵ=\frac{aM_2E}{Z_1Z_2e^2M_1}$$
(10)
and

$$a=\frac{0.885a_0}{(\sqrt{Z_1}+\sqrt{Z_2})^{2/3}},\gamma =\frac{4M_2}{M_1}$$
(11)
$`M_1=M`$ is the mass of the incident Q-ball; $`M_2`$ is the mass of the target nuclei; $`Z_1e`$ and $`Z_2e`$ are their electric charges; we assume that $`M_1>M_2`$. Nuclear losses dominates for $`\beta 10^4`$.

The energy losses of SECS in the earth mantle and earth core:

The energy losses of SECS in the earth mantle and earth core have been computed for different $`\beta `$-regions and for different charges of the Q-ball core, using the same general procedures used in the past for computing the energy losses in the earth of magnetic monopoles and nuclearites .

In Fig. 2 is presented the energy losses of Q-balls type SECS in the Earth mantle.

### 3.2 Interaction of Q-balls type SENS

The Q-ball interior of SENS is characterized by a large Vacuum Expectation Value (VEV) of certain squarks, and may be sleptons and Higgs fields . The $`SU(3)_c`$ symmetry is broken and deconfinement takes place inside the Q-ball. If a nucleon enters this region of deconfinement, it dissociates into three quarks, some of which may then become absorbed in the condensate via the reaction .

$$qq\stackrel{~}{q}\stackrel{~}{q}$$
(12)
In practice the reaction looks like

$$(Q)+Nucleon(Q+1)+pions$$
(13)
and sometimes as

$$(Q)+Nucleon(Q+1)+Kaons$$
(14)
If it is assumed that the energy released in (13) is the same as in typical hadronic processes (about 1 GeV per nucleon), this energy is carried by 2 or 3 pions (or 2 kaons). The cross section for reactions (13) and (14) is determined by the Q-ball radius $`R`$

$$\sigma 6\times 10^{34}Q^{1/2}\left(\frac{1TeV}{M_S}\right)^2cm^2$$
(15)
The corresponding mean free path $`\lambda `$ is

$$\lambda =\frac{1}{\sigma N}$$
(16)
According to References \[5-8\] the energy loss of SENS moving with velocities in the range $`10^4<\beta <10^2`$ is constant and is given by

$$\frac{dE}{dx}\frac{\zeta }{\lambda }=\zeta N10^{34}Q^{1/2}\left(\frac{1TeV}{M_S}\right)^2cm^2$$
(17)
where $`Nisthenumberofatoms/cm^3`$ and $`\zeta `$ is the energy released in the decay. The energy loss of SENS is due to the energy released from the absorbed nucleon mass. SENS lose a very small fraction of their kinetic energy and are able to traverse the Earth without attenuation for all masses of our interest.

## 4 ENERGY LOSSES IN DETECTORS

### 4.1 Light Yield of Q-balls type SECS

For SECS we distinguish two contributions to the Light Yield in scintillators: the primary Light Yield and the secondary Light Yield.

The primary Light Yield: is due to the direct excitation and ionization produced by the SECS in the medium. The energy losses in the MACRO liquid scintillator is computed from the energy losses of protons in hydrogen and carbon

$$\left(\frac{dE}{dx}\right)_{SECS}=\frac{1}{14}[2\left(\frac{dE}{dx}\right)_H+12\left(\frac{dE}{dx}\right)_C]$$
(18)
$$=SP=\frac{SL\times SH}{SL+SH}$$
(19)
where $`SP`$ is the stopping power of SECS, which reduces to $`SL`$ at low $`\beta `$ and to $`SH`$ at high $`\beta `$, so at very high $`\beta `$; the $`SP`$ stopping power coincides with the Bethe Bloch formula for electric energy losses.

1. For $`Q=1`$ the energy losses of SECS in hydrogen and carbon is computed from  adding an exponential factor due to the experimental data .

i) For $`10^5<\beta <5\times 10^3`$ we obtained the following formula

$$\left(\frac{dE}{dx}\right)_{SECS}=C[1exp(\frac{\beta }{7\times 10^4})^2]\frac{MeV}{cm}$$
(20)
where $`C=1.3\times 10^5\beta `$.
ii) For $`5\times 10^3<\beta <10^2`$ we used the following formula

$$SP=SP_H+SP_C=\left(\frac{dE}{dx}\right)_{SECS}$$
(21)
where

$$SP_H=\frac{SL_H\times SH_H}{SL_H+SH_H}$$
(22)
$$SP_C=\frac{SL_C\times SH_C}{SL_C+SH_C}$$
(23)
and

$$SL=A_1E^{0.45},SH=A_2Ln(1+\frac{A_3}{E}+A_4E)$$
(24)
where ($`A_{i=1,4}`$) are constants obtained from experimental data, and $`E`$ is the energy of a proton with velocity $`\beta `$.

2. For SECS with $`Q=Z_1e`$ the energy losses for $`10^5<\beta <10^2`$ are given by

$$\left(\frac{dE}{dx}\right)_{SECS}=F(Z_1,Z_2)[1exp(\frac{\beta }{7\times 10^2})^2]$$
(25)
where

$$F(Z_1,Z_2)=\frac{8\pi e^2a_0\beta }{\alpha }\frac{Z_1^{7/6}N_e}{(Z_1^{2/3}+Z_2^{2/3})^{3/2}}$$
(26)
where $`Z_2`$ is the atomic number of the target atom, $`N_e`$ the density of electrons and $`\alpha `$ is the fine structure constant.

The primary Light Yield of SECS is given by

$$\left(\frac{dL}{dx}\right)_{SECS}=A[\frac{1}{1+AB\frac{dE}{dx}}]\frac{dE}{dx}$$
(27)
where $`\frac{dE}{dx}`$ is the energy loss of SECS; $`A`$ and $`B`$ are parameters depending only on the velocity of SECS.

The secondary Light Yield: we considered the elastic or quasielastic recoil of hydrogen and carbon nuclei. The light yield $`L_p`$ from a hydrogen or carbon nucleus of given initial energy $`E`$ is computed as

$$L_p(E)=_0^E\frac{dL}{dx}(ϵ)S_{\text{tot}}^1𝑑ϵ$$
(28)
where $`S_{\text{tot}}`$ is the sum of electronic and nuclear energy losses. The nuclear energy losses are given in ref. . The secondary light yield is then

$$\left(\frac{dL}{dx}\right)_{\text{secondary}}=N_0^{T_m}L_p(T)\frac{d\sigma }{dT}𝑑T$$
(29)
where $`T_m`$ is the maximum energy transferred and $`\frac{d\sigma }{dT}`$ is the differential scattering cross section, given in ref. .

In Fig. 3 is presented the light yield of SECS in MACRO liquid scintillator as function of the SECS velociy $`\beta `$.

### 4.2 Energy losses of Q-balls type SECS in streamer tubes

The composition of the gas in the MACRO limited streamer tubes is 73% helium, 27% n-pentane in volume . The pressure is about one atmosphere and the resulting density is low (in comparison with the density of the other detectors): $`\rho _{gas}=0.856\text{mg/cm}^3`$. The energy losses of MMs in the streamer tubes have been discussed in .

The ionization energy losses of SECS in the streamer tubes of the MACRO experiment for $`10^3<\beta <10^2`$ was computed with the same procedure used for scintillators, but using the density and the chemical composition of streamer tubes.

For $`Q=13e`$ the energy losses are calculated from ref.  but we have omitted the exponential factor which takes into account the energy gap in organic scintillators.

The Drell effect does not occurs because SECS are not magnetically charged.

The threshold for ionizing n-pentane occur at $`\beta 10^3`$.

### 4.3 Restricted Energy Losses of SECS in the Nuclear Track Detector CR39

The quantity of interest for the CR39 nuclear track detector is the Restricted Energy Loss (REL), that is, the energy deposited within $``$ 100 Å diameter from the track.

The REL in CR39 has already been computed for MMs of $`g=g_D`$ and $`g=3g_D`$ and for dyons with $`q=e`$, $`g=g_D`$ . We have checked these calculations and extended them to other cases of interest .

The chemical composition of CR39 is $`\text{C}_{12}\text{H}_{18}\text{O}_7`$ , and the density is 1.31 $`\text{g/cm}^3`$. For the computation of the REL only energy transfers to atoms above $`12`$ eV are considered, because it is estimated that $`12`$ eV are necessary to break the molecular bonds .

At low velocities ($`3\times 10^5<\beta <10^2`$) there are two contributions to REL: the ionization and the atomic recoil contributions.

The ionization contribution was computed with Ziegler’s fit to the experimental data .

The atomic recoil contribution to REL was calculated using the interaction potential between an atom and a SECS which is equal to the electric potential  given by

$$V(r)=\frac{Z_1Z_2e^2}{r}\varphi (r)$$
(30)
where $`r`$ is the distance between the core of SECS and the target atom, $`e`$ is the electric charge of the core, $`Z_2`$ is the atomic number of the target atom. The function $`\varphi (r)`$ is the screening function given by

$$\varphi (r)=\underset{1}{\overset{3}{}}C_iexp[\frac{b_ir}{a}]$$
(31)
where $`a`$ is the screening length

$$a=0.8853\frac{a_0}{(Z_1^{\frac{1}{2}}+Z_2^{\frac{1}{2}})^{\frac{2}{3}}}$$
(32)
where $`a_0`$ is the Bohr radius; the coefficients are restricted such that

$$\underset{1}{\overset{3}{}}C_i=1$$
(33)
Assuming the validity of this potential, we calculated the relation between the scattering angle $`\theta `$ and the impact parameter b. From this relation, the differential cross section $`\sigma (\theta )`$ is readily obtained as

$$\sigma (\theta )=(db/d\theta )b/\mathrm{sin}\theta $$
(34)
The relation between the transferred kinetic energy K and the scattering angle $`\theta `$ is given by the relation

$$K=4E_{\text{inc}}\mathrm{sin}^2(\theta /2)$$
(35)
where $`E_{\text{inc}}`$ is the energy of the atom relative to the SECS in the center of mass system. The Restricted Energy Losses are finally obtained by integrating the transferred energies as

$$\frac{dE}{dx}=N\sigma (K)𝑑K$$
(36)
where $`N`$ is the number density of atoms in the medium, $`\sigma (K)`$ is the differential cross section as function of the transferred kinetic energy K.

In Fig. 4 is presented the restricted energy losses of the Q-ball type SECS in the nuclear solide track detector CR39.

## 5 CONCLUSION

Supersymmetric generalizations of the Standard Model allow for stable non-topological solitons of Q-ball type which may be considered bags of squarks and sleptons and thus have non-zero baryon and lepton numbers, as well as the electric charge \[1-3\]. These solitons can be produced in the Early Universe, can affect the nucleosynthesis of light elements, and can lead to a variety of other cosmological consequences.

In this paper, we computed the energy losses of Q-balls of type SENS and SECS. Using these energy losses and a rough model of the earth’s composition and density profiles, we have computed the geometrical acceptance of the MACRO detector for Q-balls type SECS with $`v=250km/s`$ as function of the Q-ball mass $`M`$. We have calculated the accessible region in the plane (mass, velocity) of SECS reaching MACRO from above and below.

We also presented a systematic analysis of the energy deposited in scintillators, streamer tubes and $`CR39`$ nuclear track detectors by SECS in forms useful for their detection. In particular we computed the light yield in scintillators, the ionization in streamer tubes and the REL in nuclear track detectors.

MACRO is sensitive to both SECS and SENS. A good upper flux limit may be obtained at the level of few times $`10^{16}cm^2s^1sr^1`$.

Acknowledgements : I gratefully acknowledge Prof. G. Gaicomelli for his continous disponibility and for very useful critical discussions. This work was supported by ICTP and INFN grants.


text:
# Number of distinct sites visited by 𝑁 random walkers on a Euclidean lattice

## I INTRODUCTION

Usually, the extremely successful theory of random walks is only concerned with problems that involve a single ($`N=1`$) random walker. A solid reason for this is the understanding that the average properties of the single diffusing walker serve to describe the global properties of systems formed by many walkers. However, there are other interesting diffuson problems that involve many random walkers for which the diffusive behavior of every walker of the total of $`N`$ is relevant, i.e., diffusion process that cannot be described by averaging over the properties of a single walker . The problem of evaluating the time spent by the first $`j`$ particles out of a total of $`N`$ to escape from a given region is a clear example . Another important example, which is the subject of this paper, is the problem of evaluating the average number $`S_N(t)`$ of distinct sites visited (or territory explored) by a set of $`N`$ independently diffusing random walkers up to time $`t`$ .

The case $`N=1`$ has been studied in detail since it was posed by Dvoretzky and Erdös  and is discussed in many general references . However, the multiparticle ($`N>1`$) version of this problem has been systematically treated only after the pioneering works of Larralde et al. . These authors addressed the problem of evaluating the territory covered by a set of $`N`$ independent random walkers, all initially placed at the same point, that diffuse with steps of finite variance on Euclidean lattices. They found asymptotic expressions for $`S_N(t)`$ for $`N1`$, and described the existence of three time regimes. Their results can be summarized as follows:

$$S_N(t)\{\begin{array}{cc}t^d\hfill & \hfill t<tt_\times \\ \multicolumn{2}{c}{}\\ t^{d/2}\mathrm{ln}^{d/2}\left(x\right),\hfill & \hfill t_\times tt_\times ^{^{}}\\ \multicolumn{2}{c}{}\\ NS_1(t),\hfill & \hfill t_\times ^{^{}}t\end{array},$$
(1)
where $`x=N`$ for $`d=1`$, $`x=N/\mathrm{ln}t`$ for $`d=2`$ and $`x=N/\sqrt{t}`$ for $`d=3`$ . The properties of $`S_1(t)`$ are well known; in particular, $`S_1(t)t^{1/2}`$ for $`d=1`$, $`S_1(t)t/\mathrm{ln}t`$ for $`d=2`$ and $`S_1(t)t`$ for $`d=3`$. In the very short-time regime ($`tt_\times `$), or regime I, there are so many particles at every site that all the nearest neighbors of the already visited sites are reached at the next step, so that the number of distinct sites visited grows as the volume of an hypersphere of radius $`t`$, $`S_N(t)t^d`$. The regime III ($`t_\times ^{^{}}t`$), or long-time regime, corresponds to the final stage in which the walkers move far away from each other so that their trails (almost) never overlap and $`S_N(t)NS_1(t)`$. The crossover time from regime I to regime II is given by $`t_\times \mathrm{ln}N`$ for every lattice. This can be easily understood if we take into account that the number of particles on the outer visited sites for very short times will decrease as $`N/z^t`$, where $`z`$ is the coordination number of the lattice, so that the overlapping regime will break approximately when $`N/z^t1`$ or, equivalently, $`t_\times \mathrm{ln}N`$. Regime III never appears in the one-dimensional case (i.e., $`t_\times ^{^{}}\mathrm{}`$), but $`t_\times ^{^{}}e^N`$ for $`d=2`$ and $`t_\times ^{^{}}=N^2`$ for $`d=3`$. These crossover times will be obtained readily from the mathematical formalism discussed in the present paper. The most interesting regime is regime II ($`t_\times tt_\times ^{^{}}`$), or the intermediate regime. For this time regime we will obtain explictly the main term and the first two corrective terms of the asymptotic expression of $`S_N(t)`$ for $`N1`$. Higher corrective terms could be calculated as our method allows them to be obtained in a systematic way. The contribution of these corrective terms cannot be ignored even for very large values of $`N`$ because they decay logarithmically with $`N`$. However, as we will see in section V, the use of two corrective terms leads to a very good agreement with simulation results for relatively small values of $`N`$ ($`N100`$).

The paper is organized as follows. The asymptotic evaluation of $`S_N(t)`$ for a $`d`$-dimensional Euclidean lattice is discussed in detail in Sec. II. Some geometric implications of this result are discussed in Sec. III. In Sec. IV, we compare our zeroth- (i.e., main), first- and second-order term approximation for $`S_N(t)`$ with other approximations and with computer simulations for one-, two- and three-dimensional simple Euclidean cubic lattices. The paper ends with some remarks on the applicability of this method to other diffusion problems and different media. Some technical details are discussed in an Appendix.

## II THE NUMBER OF DISTINCT SITES VISITED

We consider a group of $`N`$ random walkers starting from an origin site $`𝐫=0`$ at time $`t=0`$. A survival probability, $`\mathrm{\Gamma }_N(t,𝐫)`$, is defined as the probability that site $`𝐫`$ has not been visited by the random walkers before time $`t`$. Similarly, we can define a mortality function, $`1\mathrm{\Gamma }_N(t,𝐫)`$, as the probability that site $`𝐫`$ has been visited by at least one walker in the time interval $`(0,t)`$. The relationship between the number of distinct sites visited, $`S_N(t)`$, and the survival probability is

$$S_N(t)=\underset{𝐫}{}\left\{1\mathrm{\Gamma }_N(t,𝐫)\right\}.$$
(2)
For independent random walkers, we have $`\mathrm{\Gamma }_N(t,𝐫)=\left[\mathrm{\Gamma }_t(𝐫)\right]^N`$, where $`\mathrm{\Gamma }_t(𝐫)\mathrm{\Gamma }_1(t,𝐫)`$ is the one-particle survival probability. Next, the discrete analysis implicit in Eq. (2) is replaced by a continuous one. Thus, we write

$$S_N(t)=_0^{\mathrm{}}\left[1\mathrm{\Gamma }_t^N(r)\right]𝑑v_0r^{d1}𝑑r,$$
(3)
where $`v_0`$ is the volume (i.e., the number of lattice sites) of the hyphersphere with unit radius. It has been found for Euclidean lattices that

$$\mathrm{\Gamma }_t(r)\stackrel{~}{\mathrm{\Gamma }}_t(\xi )=1A\xi ^{2\mu }e^{d\xi ^2/2}\left(1+\underset{n=1}{\overset{\mathrm{}}{}}h_n\xi ^{2n}\right),$$
(4)
for $`\xi r/\sqrt{2dDt}1`$. Here $`D`$ is the diffusion coefficient defined through the Einstein relation $`r^22dDt`$, $`t\mathrm{}`$, with $`r^2`$ being the mean-square displacement of a single random walker. The values of $`A`$, $`\mu `$ and $`h_1`$ for $`d=1`$, $`2`$ and $`3`$ are shown in Table I. A change to the new variable $`\xi `$ and integration by parts (taking into account that $`\stackrel{~}{\mathrm{\Gamma }}_t(\mathrm{})=1`$), yields

$$S_N(t)=v_0(2dDt)^{d/2}J_N(d;0,\mathrm{}),$$
(5)
where

$$J_N(d;a,b)=_a^bN\mathrm{\Gamma }_t^{N1}(\xi )\frac{d\mathrm{\Gamma }_t(\xi )}{d\xi }\xi ^d𝑑\xi ,$$
(6)
In order to evaluate the asymptotic behavior of $`J_N(d;0,\mathrm{})`$ it is convenient to make the decomposition

$$J_N(d;0,\mathrm{})=J_N(d;0,\xi _\times )+J_N(d;\xi _\times ,\mathrm{}),$$
(7)
where $`\xi _\times `$ is a value that should satisfy the following two conditions: (a) $`\xi _\times `$ is large enough so that $`\mathrm{\Gamma }_t(r)`$ can be well approximated by its asymptotic approximation $`\stackrel{~}{\mathrm{\Gamma }}_t(\xi )`$ for $`\xi \xi _\times `$, and (b) small enough so that

$$\stackrel{~}{\mathrm{\Gamma }}_t^N(\xi _\times )=1/N^p,$$
(8)
with $`p>1`$ (say $`p=2`$). From Eq. (4) it is straightforward to see that

$$\xi _\times ^2\mathrm{ln}N$$
(9)
satisfies both conditions. On the other hand, because at most $`d\mathrm{\Gamma }/d\xi =𝒪(1)`$, and $`\mathrm{\Gamma }_t(\xi )`$ is a monotonous growing function, $`J_N(d;0,\xi _\times )`$ is bounded by a term that goes as $`N\stackrel{~}{\mathrm{\Gamma }}_t^{N1}(\xi _\times )\xi _\times ^d`$, or equivalently, from Eq. (8), by a term that goes mainly as $`N^{1p}`$. But shortly we will show that $`J_N(d;\xi _\times ,\mathrm{})`$ goes essentially as $`\mathrm{ln}^{d/2}N`$; this means that $`J_N(d;0,\xi _\times )`$ is asymptotically smaller than any term in the asymptotic expansion for $`J_N(d;0,\mathrm{})`$ and thus we can write

$$J_N(d;0,\mathrm{})J_N(d;\xi _\times ,\mathrm{}),N1.$$
(10)
The previous discussion is illustrated in Fig. 1 for the one-dimensional case. In this figure we have represented the integrand of $`J_N(1;0,\mathrm{})`$ for increasing values of $`N`$ using as survival probability the exact value $`\mathrm{\Gamma }_t(\xi )=\text{erf}(\xi /\sqrt{2})`$  and the asymptotic expression of Eq. (4) up to first order ($`n=1`$). Notice that the area below the solid \[broken\] curve is just the exact \[asymptotic approximate\] value of $`S_N(t)/(8Dt)^{1/2}=J_N(1;0,\mathrm{})`$. The value of $`\xi _\times `$ as given by Eq. (8) with $`p=2`$ is marked with a symbol. It is clear from the figure that, for large $`N`$, (a) the integrand of $`J_N(1;\xi _\times \mathrm{})`$ is well represented by the asymptotic expression of $`\mathrm{\Gamma }_t(\xi )`$, and (b) that, as stated for the general case in Eq. (10), $`J_N(1;0,\xi _\times )J_N(1;\xi _\times \mathrm{})J_N(1;0,\mathrm{})`$.

From Eq. (4), one easily finds that

$$\frac{d\stackrel{~}{\mathrm{\Gamma }}_t(\xi )}{d\xi }\left[1\stackrel{~}{\mathrm{\Gamma }}_t(\xi )\right]^1=2\xi \underset{n=0}{\overset{\mathrm{}}{}}j_n\xi ^{2n},$$
(11)
with $`j_0=d/2`$, $`j_1=\mu `$, $`j_2=h_1`$, …By inserting Eq. (11) into Eq. (6) one has the following expansion for $`J_N(d;\xi _\times ,\mathrm{})`$:

$$J_N(d;\xi _\times ,\mathrm{})2N\underset{n=0}{\overset{\mathrm{}}{}}j_nK_{N1}\left(d2n+1\right),$$
(12)
with

$$K_N(\alpha )=_{\xi _\times }^{\mathrm{}}\xi ^\alpha \stackrel{~}{\mathrm{\Gamma }}_t^N(\xi )\left[1\stackrel{~}{\mathrm{\Gamma }}_t(\xi )\right]𝑑\xi .$$
(13)
By means of the substitution

$$\stackrel{~}{\mathrm{\Gamma }}_t(\xi )=e^z,$$
(14)
we get a more convenient expression for $`K_N(\alpha )`$:

$$K_N(\alpha )=_0^{z_\times }e^{Nz}\left(e^z1\right)\xi ^\alpha \frac{d\xi }{dz}𝑑z,$$
(15)
where, from Eq. (8), $`z_x\mathrm{ln}N/N`$. The integral in (15) is of Laplace type but it is not possible to use Watson’s lemma directly to get its asymptotic behavior because $`\xi ^\alpha (d\xi /dz)`$ has a logarithmic singularity at $`z=0`$ . The evaluation of $`K_N(\alpha )`$ requires the inversion of (14) to obtain $`\xi (z)`$. By using Eqs. (4) and (14) we get

$$\frac{d}{2}\xi ^2+\mathrm{ln}A+\mu \mathrm{ln}\xi ^2+\mathrm{ln}\left(1+\underset{n=1}{\overset{\mathrm{}}{}}h_n\xi ^{2n}\right)=\mathrm{ln}\left(1e^z\right).$$
(16)
The function $`\xi (z)`$ can be readily obtained from this equation to first approximation: Notice that, as long as

$$\xi ^2|\mathrm{ln}A|,$$
(17)
the left hand side of Eq. (16) can be approximated by $`d\xi ^2/2`$, so that the first-order solution to Eq. (16) is $`\xi ^2(z)2\mathrm{ln}[1\mathrm{exp}(z)]/d`$. Equation (16) can be systematically solved in order to get higher-order approximations (see Appendix). The result is

$$\xi =x^{1/2}\underset{n=1}{\overset{\mathrm{}}{}}\delta _nx^n,$$
(18)
where $`x=(d/2)/\mathrm{ln}[1\text{exp}(z)]`$. The substitution of Eq. (18) into Eq. (15) (see Eq. (47) in the Appendix) yields

$$K_N(\alpha )=\underset{n=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{n}{}}\frac{2^{\left(\alpha 1\right)/2}}{d^{\left(\alpha +1\right)/2}}k_m^{(n)}I(\frac{\alpha }{2}n\frac{1}{2},m;N),$$
(19)
where

$$I(n,m;N)_0^{z_\times }𝑑ze^{Nz}(\mathrm{ln}z)^n\mathrm{ln}^m(\mathrm{ln}z).$$
(20)
The evaluation of $`I_N(n,m;N)`$ for $`N\mathrm{}`$ has been discussed in . For the sake of completeness, we give here explicitly their expressions up to the order required to find $`S_N(t)`$ to second order in $`1/\mathrm{ln}N`$:

$`I(n,0;N)`$ $``$ $`{\displaystyle \frac{1}{N}}\mathrm{ln}^nN\left[1+{\displaystyle \frac{n\gamma }{\mathrm{ln}N}}+{\displaystyle \frac{n(n1)}{2}}{\displaystyle \frac{\gamma ^2+\pi ^2/6}{\mathrm{ln}^2N}}+\mathrm{}\right],`$ (21)
$`I(n,1;N)`$ $``$ $`{\displaystyle \frac{1}{N}}\mathrm{ln}^nN\left[\mathrm{ln}\mathrm{ln}N\left(1+{\displaystyle \frac{n\gamma }{\mathrm{ln}N}}\right)+{\displaystyle \frac{\gamma }{\mathrm{ln}N}}+\mathrm{}\right],`$ (22)
$`I(n,2;N)`$ $``$ $`{\displaystyle \frac{1}{N}}\left(\mathrm{ln}^nN\right)\mathrm{ln}^2\mathrm{ln}N+\mathrm{}`$ (23)
where $`\gamma 0.577215`$ is the Euler constant. Using these results we get from Eqs. (5), (12) and (19) the following expansion for the average number of distinct sites visited on a Euclidean lattice of dimension $`d`$:

$$S_N(t)\widehat{S}_N(t)(1\mathrm{\Delta })$$
(24)
with

$`\widehat{S}_N(t)`$ $`=`$ $`v_0\left(4Dt\mathrm{ln}N\right)^{d/2},`$ (25)
$`\mathrm{\Delta }\mathrm{\Delta }(N,t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{ln}^nN{\displaystyle \underset{m=0}{\overset{n}{}}}s_m^{(n)}\mathrm{ln}^m\mathrm{ln}N`$ (26)
and where, up to second order ($`n=2`$),

$`s_0^{(1)}`$ $`=`$ $`d\omega ,`$ (27)
$`s_1^{(1)}`$ $`=`$ $`d\mu ,`$ (28)
$`s_0^{(2)}`$ $`=`$ $`d\left(1{\displaystyle \frac{d}{2}}\right)\left({\displaystyle \frac{\pi ^2}{12}}+{\displaystyle \frac{\omega ^2}{2}}\right)d\left({\displaystyle \frac{dh_1}{2}}\mu \omega \right),`$ (29)
$`s_1^{(2)}`$ $`=`$ $`d\left(1{\displaystyle \frac{d}{2}}\right)\mu \omega d\mu ^2,`$ (30)
$`s_2^{(2)}`$ $`=`$ $`{\displaystyle \frac{d}{2}}\left(1{\displaystyle \frac{d}{2}}\right)\mu ^2.`$ (31)
Here $`\omega =\gamma +\mathrm{ln}A+\mu \mathrm{ln}(d/2)`$, and $`A`$, $`\mu `$ and $`h_1`$ are given in Table I for $`d=1`$, $`2`$ and $`3`$. Notice that the time dependence of $`\mathrm{\Delta }(N,t)`$ comes from the term $`\omega `$ through the function $`A(t)`$. However, this function does not depend on time for the one-dimensional case and thus $`\mathrm{\Delta }`$ only depends on $`N`$.

Recently, Sastry and Agmon  have found an approximate formula for $`S_N(t)`$ for the one-dimensional case. The straightforward method used by these authors is based on the fact that the function $`\mathrm{\Gamma }_t^N(r)`$ that appears in the integrand of Eq. (2) approaches a step function when $`N\mathrm{}`$. In this way they found

$`S_N(t)`$ $``$ $`4\sqrt{Dt}\sqrt{\mathrm{ln}N\mathrm{ln}\sqrt{\alpha \mathrm{ln}N}}`$ (32)
$``$ $`4\sqrt{Dt\mathrm{ln}N}\left[1{\displaystyle \frac{1}{4}}{\displaystyle \frac{\mathrm{ln}\mathrm{ln}N\mathrm{ln}\alpha }{\mathrm{ln}N}}+𝒪\left({\displaystyle \frac{\mathrm{ln}^2\mathrm{ln}N}{\mathrm{ln}^2N}}\right)\right]`$ (33)
where $`\alpha `$ is given by $`\alpha =\pi \mathrm{exp}(2/\pi )1.66`$. It is instructive to compare this formula with the first-order approximation of Eq. (24) for the one-dimensional case:

$$S_N(t)4\sqrt{Dt\mathrm{ln}N}\left[1\frac{1}{4}\frac{\mathrm{ln}\mathrm{ln}N2\omega }{\mathrm{ln}N}+𝒪\left(\frac{\mathrm{ln}^2\mathrm{ln}N}{\mathrm{ln}^2N}\right)\right],$$
(34)
Note that the prefactor $`4(Dt\mathrm{ln}N)^{1/2}`$ of the formula of Sastry and Agmon is in agreement with that of Eq. (34). In Ref. , they found it “amusing” that the value $`\alpha =1`$ produces very good agreement between the approximation of Eq. (32) and the exact numerical integration. Our Eq. (34) enlightens this point: Comparing Eqs. (33) and Eq. (34) one sees that $`\mathrm{ln}\alpha `$ is playing the role of $`2\omega `$. But $`\omega =\gamma \frac{1}{2}\mathrm{ln}\pi =0.0048507\mathrm{}`$ for the one-dimensional lattice, so that $`\mathrm{ln}\alpha `$ when $`\alpha =1`$ leads to a good approximation to the rigorous coefficient $`2\omega `$. The equation of Sastry and Agmon for $`\alpha =1`$ and our first-order approximation should thus be very close. This is clearly confirmed in Fig. 4.

A question to be answered is why Eq. (24) is valid for time regime II only, i.e., why it is not always valid for arbitrarily large values of time. The reason is that our formulas have been obtained by assuming that the condition (17) holds for those values of $`\xi `$ which are inside the integration interval $`[\xi _\times ,\mathrm{}]`$ of the relevant integral $`J_N(d;\xi _\times ,\mathrm{})`$ that is responsible for the asymptotic behavior of $`S_N(t)`$. This implies that for our procedure to work, it is necessary that $`\xi _\times ^2\mathrm{ln}A`$ or, from Eq. (9), that

$$\mathrm{ln}N|\mathrm{ln}A|.$$
(35)
Thus we can estimate the time $`\tau _\times `$ for which our method breaks down by solving $`|\mathrm{ln}A(\tau _\times )|\mathrm{ln}N`$. From the expressions for $`A`$ quoted in Table I one finds that $`\tau _\times e^N`$ for $`d=2`$ and $`\tau _\times N^2`$ for $`d=3`$. For $`d=1`$ and large $`N`$, the condition (35) always holds because $`A=(2/\pi )^{1/2}`$ is a constant and then $`\tau _\times =\mathrm{}`$. We see that the upper times $`\tau _\times `$ beyond which Eq. (24) is no longer valid coincide with the crossover times $`t_\times ^{^{}}`$ defined in Sec. I, so we can say that the expressions for $`S_N(t)`$ given in this paper are valid only in the time regime II. This means that our procedure marks its own limit of validity as that of regime II, thus predicting the existence of a crossover time in a natural way, i.e., as a consequence of the mathematical formalism.

## III GEOMETRIC PROPERTIES OF THE EXPLORED REGION

In this section we will give a geometric interpretation of the main result of this paper, namely, Eq. (24). The quantity $`S_N(t)`$ is by definition the volume of the region $`\mathrm{\Omega }`$ explored by $`N`$ random walkers after a time $`t`$ from their initial deposition on a given site of the lattice (if the length of the lattice bonds is taken as the unit). For very short times (regime I or $`t\mathrm{ln}N`$) the exploration is performed in a compact way because all the neighbor sites of any visited site are always visited at the next time step. Therefore, the explored region $`\mathrm{\Omega }`$ is an hypersphere whose radius grows ballistically and its volume is proportional to $`t^d`$. After the regime II is reached, the development of two qualitatively different zones in the explored volume is observed: (i) an hyperspherical compact core of visited sites, and (ii) a corona of dendritic nature characterized by filaments created by those relatively few walkers that are wandering in the outer regions, i.e., wandering at distances significantly larger than the root-mean-square displacement $`r^2^{1/2}=\sqrt{2dDt}`$ of a single walker. Figure 2 shows a snapshot of the set of sites visited by $`N=1000`$ random walkers at time $`t=900`$ (every walker makes a jump at each time unit) for dimension two. The visited sites are in white and the inner black and outer white circles delimit the corona. The radius $`R_+`$ of the outer circle is equal to the maximum displacement from the origin reached by any of the walkers at time $`t`$. It has been argued in  that the volume of this outer circle is on average given by the main term of (24), i.e., by $`\widehat{S}_N(t)=v_0(4Dt\mathrm{ln}N)^{d/2}`$. From this statement we can draw two conclusions: First, that the average radius of this outer circle is

$$R_+\left(4Dt\mathrm{ln}N\right)^{1/2},$$
(36)
and second, that the asymptotic corrective terms (given by $`\mathrm{\Delta }`$) to $`S_N(t)`$ account for the number of unvisited sites that are inside the corona. In other words, $`\mathrm{\Delta }`$ is the fraction of the volume inside the external circumference that has not been visited by any of the $`N`$ random walkers. This result can be used  to easily estimate that the thickness of the dendritic corona is approximately given by $`R_+\mathrm{\Delta }`$.

It is also noteworthy that $`\mathrm{\Delta }(N,t)`$ depends on $`t`$ very smoothly in the time regime II as this dependence is due to terms proportional to powers of $`\mathrm{ln}A(t)`$ \[and $`A(t)`$ does not change exponentially: see Table I\]. For the two-dimensional case, this statement is especially valid because $`A(t)1/\mathrm{ln}t`$. Therefore, the ratio (given by $`\mathrm{\Delta }`$) between the radial size of the corona of $`\mathrm{\Omega }`$ and the radial size of $`\mathrm{\Omega }`$ itself remains almost constant throughout the time regime II. This implies that a conveniently scaled sequence of snapshots of the set of visited sites should be very similar (in a statistical sense), i.e., we find that $`\mathrm{\Omega }`$ grows, to a large extent, in a self-similar way inside time regime II. This property is illustrated in Fig. 3. As Eq. (36) shows, the appropriate scale factor must be proportional to $`\sqrt{t}`$. This “almost” self-similar behavior disappears as the regime III is approached because the correction to the main term of $`S_N(t)`$ becomes as large as this main term, i.e., because $`\mathrm{\Delta }(N,t)`$ approaches the value 1. This transition takes place when $`t\tau _\times `$ as follows from (26), i.e., this value coincides with the threshold for regime III deduced in the previous section. From the geometric point of view this transition corresponds to the breaking of the self-similar growing behavior by the appearance of a corona of filaments as large as the compact core, which finally gives rise to a set of separated trails that (almost) never more overlap. For the two-dimensional case the transition time from regime II to regime III is so great for any significant number of walkers that it can not be studied by numerical simulation.

## IV NUMERICAL RESULTS

We carried out numerical simulations for the number of distinct sites visited by $`N=2^m`$, with $`m=0,1,\mathrm{},14`$ in two and three dimensions. For the one-dimensional case it is not necessary to carry out simulations because the survival probability is exactly known on this lattice, $`\mathrm{\Gamma }_t(𝐫)=\text{erf}\left(\xi /\sqrt{2}\right),`$ and therefore the integral for $`S_N(t)`$ as given by Eq. (3) can be computed numerically.

In our simulations, the random walkers are placed initially at the center of an hypercubic box of side $`L`$. The regime II is reached almost immediately with the number of random walkers we have used ($`t_\times 10`$ for $`N=2^{14}`$). The simulations were carried out only to a maximum time $`t=200`$ which is sufficient for the stabilization of regime II conditions. The square box side for $`d=2`$ was taken to be $`L=400`$ to avoid any random walker reaching the edge of the box before the maximum time $`t=200`$. Memory limitations forced us to reduce the box side to $`L=200`$ for the three-dimensional case. While this implies a possible appearance of finite-size effects, we can consider them as to be negligible because the average displacement of the random walkers at the maximum time is small compared with $`L/2`$. Each experiment was repeated $`10^4`$ times in order to achieve reasonable statistics.

Results are plotted in Fig. 4 for one, two and three dimensions. The dots are the simulation results (numerical integration results in one dimension) and the broken and solid lines are the prediction of Eq. (24) to first and second order, respectively. The crosses are the results of Sastry and Agmon  given by Eq. (32) with $`\alpha =1`$. The dotted lines correspond to the result of Larralde et al. given by Eq. (1) using the correct amplitude of the main term (see ). The quantity plotted is

$$𝐒\frac{1}{d}\left[\frac{S_N}{\widehat{S}_N}\right]^{2/d},$$
(37)
versus $`1/\mathrm{ln}N`$. From Eq. (24) one sees that the theoretical prediction for this quantity is $`𝐒(1/d)(1\mathrm{\Delta })^{2/d}`$. The agreement between the second-order approximation and the simulations is found to be excellent for $`N100`$. Good agreement for lower values of $`N`$ would be expected if higher-order terms in the series were included. The importance of the corrective terms is evident. For example, for the one-dimensional case, we would need to use values of $`N`$ as large as $`10^{25}`$ in order to obtain the same precision with the main term as we get with the main and two corrective terms for values of $`N`$ as small as $`2^6`$. Similar statements can be made for the other lattices, as Fig. 4 shows.

## V REMARKS

In this paper we have developed a method for calculating the mean number of distinct sites visited by $`N`$ independent random walkers on Euclidean lattices. The method allows the systematic calculation of the main and corrective asymptotic terms to any order for large $`N`$. These corrective terms are generally non-negligible as they (essentially) decay as powers of $`1/\mathrm{ln}N`$. However, we found that the main and first two corrective terms lead to reasonably good results when relatively small values of $`N`$ are used (say, for $`N2^7`$). In Sec. III we proposed a geometric meaning for the main and corrective terms: the main term would account for the volume of the set of visited sites if the exploration of the random walkers were compact, and the corrective terms just improve this rough estimate because, in the outer regions, the exploration performed by the (relatively few) random walkers that move there is really not compact, thus leading to the formation of a non-compact (a dendritic) external ring in the set of visited sites. We hope the above results and ideas could serve as a basis to gain insight into problems with interacting random walkers.

The method developed here for calculating $`S_N(t)`$ is also useful for evaluating other statistical quantities related to the diffusion of a set of independent random walkers. An example is the number $`S_{N+}(t)`$ of sites visited by $`N`$ random walkers on an one-dimensional lattice along a given direction . It turns out that the moments (of arbitrary order) of $`S_{N+}(t)`$ can be readily obtained through a slight modification of Eq. (24). Another example is the first passage time $`t_{1,N}(r)`$ to a distance $`r`$ of the first random walker of a set of $`N`$. First passage times are relevant statistical quantities in the study of diffusion processes where the arrival of the first particles at a given site produces a significant effect (a “trigger” effect). These quantities have been calculated for one dimension  (and for some classes of fractals ) but little is known for dimensions greater than one . The approximate compact form of the set of visited sites allows one to estimate the first passage time via the relation $`S_N(t_{1,N}(r))v_0r^d`$ , which means geometrically that we consider the region inside the hypersphere of radius $`r`$ where a random walker has arrived by time $`t_{1,N}(r)`$ as completely visited (a compact exploration in the sense of de Gennes ). Results on $`S_{N+}(t)`$ and $`t_{1,N}(r)`$ obtained using the above ideas will be reported elsewhere.

The function $`S_N(t)`$ we have studied is indeed an important quantity concerning the diffusion of $`N`$ independent random walkers but there are still many open questions in this problem. One can think, for example, of the absorption probability of the set of $`N`$ random walkers on a lattice with a random distribution of point-like traps. This problem can be formulated in terms of the moments of the number of distinct sites visited by the set of $`N`$ walkers. A prediction for the variance of the number of visited sites is a necessary requisite to tackle this interesting problem as the first-order approximation based only on its first moment, i.e., on $`S_N(t)`$, seems to be very imprecise . As no relationship is known for moments of order higher than one, the absorption problem remains unsolved.

Finally, it should be pointed out that the expression for $`S_N(t)`$ given in this paper can be extended to fractal media with some slight changes. We are currently running simulations for deterministic (Sierpinski gasket) and stochastic (percolation aggregate) fractals. Results for these substrates will be published elsewhere.

###### Acknowledgements.

Partial support from the DGICYT (Spain) through grant no PB97-1501 and from the Junta de Extremadura-Fondo Social Europeo through grant IPR99C031 is acknowledged.

##

We will show in this Appendix how to get Eq. (19) from Eq. (15). Let us start by showing that the solution $`\xi (z)`$ of Eq. (16) for $`z0`$ has the form given in Eq. (18). For simplicity of notation we will write $`u=\xi ^2`$, $`\varphi =1\text{exp}(z)`$ and $`c=d/2`$. Hence, Eq. (16) takes the form

$$\frac{c}{u}+\mu \mathrm{ln}u+\mathrm{ln}A+\mathrm{ln}\left(1+\underset{n=1}{\overset{\mathrm{}}{}}h_nu^n\right)=\mathrm{ln}\varphi .$$
(38)
In the limit $`z0`$, it is clear that $`u0`$ and $`\varphi 0`$. This means that, as long as $`1/u|\mathrm{ln}(A)|`$, the first term on the right-hand side of (38) is the most divergent one so that, as a first approximation, we have

$$u\frac{c}{\mathrm{ln}\varphi }x.$$
(39)
This first-order approximation was already obtained in Sec. II \[see below Eq. (17)\]. A better approximation is achieved by writing $`u=x(1+ϵ)`$, with $`ϵ`$ an small quantity. The substitution of this expression in (38) yields

$$ϵϵ^2+\frac{\mu x}{c}\mathrm{ln}x+\frac{x}{c}\mathrm{ln}A+\frac{\mu x}{c}ϵ\frac{\mu x}{2c}ϵ^2+\frac{h_1x^2}{c}+\frac{h_1x^2ϵ}{c}+\mathrm{}=0,$$
(40)
where (39) has been taken into account. This equation can be solved by writing $`ϵ`$ as

$$ϵ=\underset{n=1}{\overset{\mathrm{}}{}}ϵ_nx^n,$$
(41)
and inserting it in (40). We thus find the following values for $`ϵ_n`$ up to $`n=2`$:

$`ϵ_1`$ $`=`$ $`{\displaystyle \frac{1}{c}}\mathrm{ln}\left(Ax^\mu \right),`$ (42)
$`ϵ_2`$ $`=`$ $`{\displaystyle \frac{1}{c^2}}\mathrm{ln}^2\left(Ax^\mu \right)+{\displaystyle \frac{\mu }{c^2}}\mathrm{ln}\left(Ax^\mu \right){\displaystyle \frac{h_1}{c}}.`$ (43)
Therefore

$$\xi (z)=u^{1/2}=x^{1/2}\left(1+ϵ\right)^{1/2}=x^{1/2}\underset{n=0}{\overset{\mathrm{}}{}}\delta _nx^n,$$
(44)
where $`\delta _0=1`$ and

$`\delta _1`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}\left(Ax^\mu \right)}{2c}},`$ (45)
$`\delta _2`$ $`=`$ $`{\displaystyle \frac{1}{8c^2}}\left[\mathrm{ln}^2\left(Ax^\mu \right)+4\mu \mathrm{ln}\left(Ax^\mu \right)4ch_1\right].`$ (46)
The evaluation of the integral for $`K_N(\alpha )`$ in (15) requires the expression of $`\xi ^\alpha d\xi /dz`$ as a function of $`z`$. From Eq. (44) and taking into account that $`d\xi /dz=(d\xi /dx)(dx/dz)`$ and $`dx/dz=[x^2/(c\varphi )]d\varphi /dz`$, we find that

$$\xi ^\alpha \frac{d\xi }{dz}=\frac{1}{2c\varphi }\frac{d\varphi }{dz}x^{(1\alpha )/2}\left[1+\underset{n=1}{\overset{\mathrm{}}{}}x^n\underset{m=1}{\overset{n}{}}\widehat{k}_m^{(n)}\mathrm{ln}^m\left(Ax^\mu \right)\right],$$
(47)
where the coefficients $`\widehat{k}_m^{(n)}`$, $`m=0,\mathrm{},n`$ for $`n=1`$, $`2`$ are

$`\widehat{k}_0^{(1)}`$ $`=`$ $`{\displaystyle \frac{\mu }{c}},`$ (48)
$`\widehat{k}_1^{(1)}`$ $`=`$ $`{\displaystyle \frac{\alpha 1}{2c}},`$ (49)
$`\widehat{k}_0^{(2)}`$ $`=`$ $`{\displaystyle \frac{(\alpha 3)h_1}{2c}}+{\displaystyle \frac{\mu ^2}{c^2}},`$ (50)
$`\widehat{k}_1^{(2)}`$ $`=`$ $`{\displaystyle \frac{\mu (2\alpha )}{2c^2}},`$ (51)
$`\widehat{k}_2^{(2)}`$ $`=`$ $`{\displaystyle \frac{\alpha (\alpha 4)+3}{8c^2}}.`$ (52)
Let us use $`\widehat{K}_N(\alpha ,z)`$ to denote the integrand of Eq. (15), i.e.,

$$K_N(\alpha )=_0^{z_\times }\widehat{K}_N(\alpha ,z)𝑑z.$$
(53)
Then, from Eq. (47),

$$\widehat{K}_N(\alpha ,z)=\frac{1}{2c}e^{Nz}e^zx^{(1\alpha )/2}\left[1+\underset{n=1}{\overset{\mathrm{}}{}}x^n\underset{m=1}{\overset{n}{}}\widehat{k}_m^{(n)}\mathrm{ln}^m\left(Ax^\mu \right)\right].$$
(54)
Writing $`e^z=1+𝒪(z)`$, $`x=(c/\mathrm{ln}z)[1+𝒪(z/\mathrm{ln}z)]`$ and $`\mathrm{ln}(Ax^\mu )=\mathrm{ln}A\mu \mathrm{ln}(\mathrm{ln}z)+\mu \mathrm{ln}c+𝒪(z/\mathrm{ln}z)`$, Eq. (54) becomes

$$\widehat{K}_N(\alpha ,z)=[1+𝒪(z)]\frac{1}{2c^{(\alpha +1)/2}}e^{Nz}(\mathrm{ln}z)^{(\alpha 1)/2}\underset{n=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{n}{}}k_m^{(n)}(\mathrm{ln}z)^n\mathrm{ln}^m(\mathrm{ln}z)$$
(55)
where the coefficients $`k_m^{(n)}`$ up to second order ($`n=2`$) are

$`k_0^{(1)}`$ $`=`$ $`(\alpha 1){\displaystyle \frac{\omega }{2}}\mu ,`$ (56)
$`k_1^{(1)}`$ $`=`$ $`(1\alpha ){\displaystyle \frac{\mu }{2}},`$ (57)
$`k_0^{(2)}`$ $`=`$ $`(3\alpha )(1\alpha ){\displaystyle \frac{\omega ^2}{8}}+\mu (2\alpha )\omega +\mu ^2+{\displaystyle \frac{h_1c}{2}}\left(\alpha 3\right),`$ (58)
$`k_1^{(2)}`$ $`=`$ $`\mu \left[(\alpha 2)\mu +(\alpha 3)(1\alpha ){\displaystyle \frac{\omega }{4}}\right],`$ (59)
$`k_2^{(2)}`$ $`=`$ $`{\displaystyle \frac{\mu ^2}{8}}(\alpha 3)(\alpha 1),`$ (60)
and $`\omega =\gamma +\mathrm{ln}A+\mu \mathrm{ln}c`$. Finally, inserting Eq. (55) into Eq. (53) we get Eq. (19). It should be noted that we have approximated the factor $`1+𝒪(z)`$ of Eq. (55) by 1. This can be done safely because the contribution of the neglected terms to the asymptotic behavior of $`K_N(\alpha )`$ decays as least as $`(\mathrm{ln}N)^{(\alpha 1)/2}/N^2`$, i.e., decays to zero faster than the contribution of the retained terms by (roughly) a factor $`N`$ \[see Eqs. (19)-(23)\].


text:
# ASCA and BeppoSAX observations of the peculiar X–ray source 4U1700+24/HD154791

## I Introduction

In optical astronomy the identification of a binary system comes in most cases from the observation of photometric and/or radial velocity variations. As not all X-ray binaries have known optical counterparts, a further effective criterium in galactic X–ray astronomy for the identification of a binary system with an accreting compact object was often based on the observed X–ray luminosity. For X–ray binaries harbouring a neutron star or possibly a black hole, luminosities LX of the order of 1034 – 1035 erg s-1 are easily reached. The diagnosis of the presence of a neutron star in most cases is directly confirmed by the observation of pulsations or thermonuclear bursts, apart from bright persistent Low Mass X–Ray Binaries (LMXRBs). X–ray binaries harbouring white dwarfs also show some distinctive features. As an example in polars and intermediate polars optical and UV observations often reveal the distinctive signatures of the presence of a white dwarf in the system. Orbital periods and light curves also add unambiguous and reliable evidence of the presence of white dwarfs in this class of X–ray binaries.
For a number of X–ray sources the identification of a class or even the diagnosis of binarity is rather difficult, especially when the observed X–ray luminosity is $`10^{33}`$ erg s-1 .
4U1700+24/HD154791 belongs to this class. The optical counterpart was identified by Garcia et al. garcia  as a late type giant on the basis of the positional coincidence with a HEAO1–A3 error box. The optical spectrum of this giant looks quite normal ddf ; garcia , even if Gaudenzi and Polcaro gaudenzi  find some interesting and variable features in its spectrum. Variable UV line emission was detected ddf ; garcia  in different IUE pointings, showing at last some unusual features in the emission from this otherwise normal giant. These high excitation lines are likely linked to the same mechanism that produces the observed X–ray emission. In spite of various attempts, no evidence of a binary orbit was obtained from radial velocity analysis of optical spectra.
The X–ray source shows extreme erratic variability, but no pulsations were detected. The rapid (10–1000 s) time variability is strongly suggestive of turbulent accretion, often observed in X–ray binaries. The X–ray spectrum is rather energetic and was measured up to 10 keV. The X–ray luminosity LX $`10^{33}`$ erg s-1 at an assumed distance of 730 pc garcia  may be marginally consistent with coronal emission, even if an evolved giant is not expected to be a strong X-ray emitter. Therefore the picture emerging from observations gives only hints in favour of a binary system, given that no “classical” feature to be associated to the presence of a compact object was found.
We have observed this source for $``$15 years both with X–ray satellites (EXOSAT, ASCA and BeppoSAX ) and with ground optical observations from the Loiano 1.5 m and 0.6 m telescopes of the Bologna Astronomical Observatory. Here we report on the ASCA and BeppoSAX observations, performed respectively on March 8, 1995 and on March 27, 1998. We also report on photometric optical UBVRI monitoring.

## II Observations

In Figure 1 we show the observed 1.5–9 keV count rate from the GIS2 and GIS3 instruments on board ASCA and the 1.5–10 keV observed count rate from the MECS2 and MECS3 instruments on board BeppoSAX . A clear increase of the 4U1700+24 count rate was detected in November 1997 by RXTE/ASM (http://space.mit.edu/XTE/ASM\_lc.html ). The BeppoSAX observation was performed approximately five months after this event, when the source had already recovered its quiescent flux. The substantial erratic variability already detected with EXOSAT ddf  is clearly present also in both observations. The source flux in the BeppoSAX observation is significantly lower than that in the ASCA observation.

The erratic source variability is clearly visible in the Power Spectral Density (PSD) shown in Figure 2, calculated on the time series of GIS2 and GIS3 count rate binned on 0.1 s. The spectra were calculated for runs with typical duration of 3000s. The PSD shown in Figure 2 is obtained by averaging the spectra of different runs and by summing adjacent frequencies with a logarithmic rebinning.

The observed X–ray source luminosity (2–10 keV) was $`L_X=1.7\times 10^{33}`$ erg s-1 in the ASCA observation and $`L_X=6\times 10^{32}`$ erg s-1 in the BeppoSAX observation assuming a distance of $``$700 pc garcia .
The X–ray energy spectrum cannot be fitted by simple single component models. The high energy ($`>`$2 keV) spectrum can be fitted by an absorbed thermal continuum, but the extrapolation of such a model at lower energies lies significantly below the measured spectrum, both in ASCA and in BeppoSAX observations.
For a thermal model, similar to that used in the low luminosity source $`\gamma `$ Cas (a suspected Be/white dwarf binary kubo ; owens ), the addition of a complex absorber (e.g. a partial absorber) is needed to model the low energy part of the spectrum. The lack of Fe emission line however requires a very low Fe abundance.
As an example the count rate spectra from the ASCA and BeppoSAX observations fitted with an optically thin thermal bremsstrahlung continuum with partial absorber (“bremss”  and “pcfabs”  models in XSPEC) are shown in Figure 3. The BeppoSAX spectrum is softer than that observed with ASCA , and very similar to that observed with EXOSAT ddf  at almost exactly the same flux level of the BeppoSAX observation.

Optical observations were performed at the Loiano 1.5 m telescope of the Bologna Astronomical Observatory during the last 15 years. HD154791, the optical counterpart of the X–ray source, is a M2–M3 giant garcia ; gaudenzi  with a rather normal optical spectrum. A simple comparison with M1–M3 III templates shows a close match with the M2 template of HD104216. In Figure 4 we report the long term UBVRI photometry of HD154791. No clear long term trend is visible. Some variability is present, in particular in the U measurements, that may be intrinsic to the source. The long-term spectral/photometric monitoring of the source is continuing.

## III Discussion

The observations we report still cannot be used to perform a “classical” diagnosis of binarity. We nevertheless note some interesting similarities with other low luminosity X–ray sources. In particular some interesting similarities can be found with the X–ray emission from $`\gamma `$ Cas. The power spectrum is strikingly similar and the energy spectrum shows a similar shape, even if no iron line is detected in 4U1700+24 .

However this close resemblance of the properties of the X–ray emission does not help to determine the presence of a compact object in a binary system, as for $`\gamma `$ Cas itself the diagnosis of binarity is not completely assessed. In fact Owens et al. owens  favour the hypothesis of a WD binary, but a completely different point of view is based on recent UV/X–ray observations of $`\gamma `$ Cas smith . Smith et al. support the hypothesis that the X–ray emission of $`\gamma `$ Cas comes from continuous flaring from the Be star. This hypothesis cannot be easily adapted to the case of 4U1700+24/HD154791 , as the much colder M giant star should not be expected to have strong and persistent X–ray flaring activity. If this is the case, i.e. the observed X–ray emission from 4U1700+24 is coronal, HD154791 should be an exception in its own class. If the similarity of the properties of the X–ray emission from 4U1700+24 and $`\gamma `$ Cas comes from a common origin, we suggest that the WD binary hypothesis is much more comfortable and more easily met.

Acknowledgements. This research is supported by the Agenzia Spaziale Italiana (ASI) and the Consiglio Nazionale delle Ricerche (CNR) of Italy. BeppoSAX is a joint program of ASI and of the Netherlands Agency for Aerospace Programs (NIVR). The ASCA observation was performed as part of the joint ESA/Japan scientific program. CB, AG and AP acknowledge a grant from “Progetti di ricerca ex-quota 60%” of Bologna University.


text:
# Simple denoising algorithm using wavelet transform

11footnotetext: e-mail address for correspondence: ravi@che.ncl.res.in

keywords: Noise, Discrete wavelet transform, Chaos, Differentiation

Application of wavelets and multiresolution analysis to reaction engineering systems from the point of view of process monitoring, fault detection, systems analysis etc. is an important topic and of current research interest (see, Bakshi and Stephanopoulos, 1994; Safavi et. al., 1997; Luo et. al., 1998; Carrier and Stephanopoulos, 1998). In the present paper we focus on one such important application, where we propose a new and simple algorithm for the reduction of noise from a scalar time series data. Presence of noise in a time–varying signal restricts one’s ability to obtain meaningful information from the signal. Measurement of correlation dimension can get affected by a noise level as small as 1% of signal, making estimation of invariant properties of a dynamical system, such as the dimension of the attractor and Lyapunov exponents, almost impossible (Kostelich and Yorke, 1988). Noise in experimental data can also cause misleading conclusions (Grassberger et. al., 1991). A host of literature exists on various techniques for noise reduction (Kostelich and Yorke, 1988; Härdle, 1990; Farmer and Sidorowich, 1991; Sauer, 1992; Cawley and Hsu, 1992; Cohen, 1995; Donoho and Johnstone, 1995; Kantz and Schreiber, 1997). For instance, Fast Fourier Transform (FFT) reduces noise effectively in those cases where the frequency distribution of noise is known (Kostelich and Yorke, 1988; Cohen, 1995; Kantz and Schreiber, 1997); singular value analysis methods (Cawley and Hsu, 1992) project the original time–series onto an optimal subspace, whereby noise components are left behind in the remaining orthogonal directions, etc. In the existing wavelet–based denoising methods (Donoho and Johnstone, 1995) two types of denoising are introduced: linear denoising and nonlinear denoising. In linear denoising, noise is assumed to be concentrated only on the fine scales and all the wavelet coefficients below these scales are cut off. Nonlinear denoising, on the other hand, treats noise reduction by either cutting off all coefficients below a certain threshold (so called ‘hard–thresholding’), or reducing all coefficients by this threshold (so called ‘soft–thresholding’). The threshold values are obtained by statistical calculations and has been seen to depend on the standard deviation of the noise (Nason, 1994).

The noise reduction algorithm that we propose here makes use of the wavelet transform (WT) which in many ways complements the well known Fourier Transform (FT) procedure. We apply our method, firstly, to three model flow systems, viz. Lorenz, Autocatalator, and Rössler systems, all exhibiting chaotic dynamics. The reasons for choosing these systems are the following: Firstly, all of them are simplified models of well–studied experimental systems. For instance, Lorenz is a simple realization of convective systems (Lorenz, 1963), while the Autocatalator and Rössler have their more complicated analogs in chemical multicomponent reactions (Rössler, 1976; Lynch, 1992). Secondly, chaotic dynamics is extremely nonlinear, highly sensitive, possesses only short–time correlations and is associated with a broad range of frequencies (Guckenheimer and Holmes, 1983; Strogatz, 1994). Because of these properties it is well known that FT methods are not applicable in a straightforward way to chaotic dynamical systems (Abarbanel, 1993). On the other hand, WT methods are particularly suited to handle not only nonlinear but nonstationary signals (Strang and Nguyen, 1996). This is because the properties of the data are studied at varying scales with superior time localization analysis when compared to FT technique. Our noise reduction algorithm is advantageous, because, as shall be shown, the threshold level for noise is identified automatically. In this study, we have used the discrete analog of the wavelet transform (DWT) which involves transforming a given signal with orthogonal wavelet basis functions by dilating and translating in discrete steps (Daubechies, 1990; Holschneider, 1995). For study purposes we corrupt one variable $`x(t)`$ for each of these systems with noise of zero mean, and then apply our algorithm for denoising. We analyze the performance of this method in all the three systems for a wide range of noise strengths, and show its effectiveness. Importantly, we then validate the applicability of the method to experimental data obtained from two chemical systems. In one system the time series data was obtained from pressure fluctuation measurements of the hydrodynamics in a fluidized bed. In the other the conductivity measurements in a liquid surfactant manufacturing experiment were analyzed.

Methodology

The noise reduction algorithm based on DWT consists of the following five steps:

Step 1:

In first step, we differentiate the noisy signal $`x(t)`$ to obtain the data $`x_d(t)`$, using the method of central finite differences with fourth order correction to minimize error (Constantinides, 1987), i.e.,

$`x_d(t)={\displaystyle \frac{dx(t)}{dt}}.`$ (1)
Step 2:

We then take DWT of the data $`x_d(t)`$ and obtain wavelet coefficients $`W_{j,k}`$ at various dyadic scales $`j`$ and displacements $`k`$. A dyadic scale is the scale whose numerical magnitude is equal to 2 (two) raised to an integer exponent, and is labeled by the exponent. Thus, the dyadic scale $`j`$ refers to a scale of magnitude $`2^j`$. In other words, it indicates a resolution of $`2^j`$ data points. Thus a low value of $`j`$ implies finer resolution while high $`j`$ analyzes the signal at coarser resolution. This transform is the discrete analog of continuous WT (Holschneider, 1995), and is given by the formula

$`W_{j,k}={\displaystyle _{\mathrm{}}^+\mathrm{}}x_d(t)\psi _{j,k}(t)𝑑t,`$ (2)
with

$`\psi _{j,k}(t)=\mathrm{\hspace{0.33em}\hspace{0.33em}2}^{j/2}\psi (2^jtk)`$

where $`j,k`$ are integers. As for the wavelet function $`\psi (t)`$ we have chosen Daubechies compactly supported orthogonal function with four filter coefficients (Daubechies, 1990; Press et. al., 1996).

Step 3:

In this step we estimate the power $`P_j`$ contained in different dyadic scales $`j`$, via

$`P_j(x)={\displaystyle \underset{k=\mathrm{}}{\overset{+\mathrm{}}{}}}|W_{j,k}|^2(j=1,2,\mathrm{})`$ (3)
By plotting the variation of $`P_j`$ with $`j`$, we see that it is possible to identify a scale $`j_m`$ at which the power due to noise falls off rapidly. This is important because as we shall see from the studies of the case examples that it provides a mean for automation in detection of threshold. The identification of the scale $`j_m`$ at which power due to noise shows the first minimum allows us to reset all $`W_{j,k}`$ upto scale $`j_m`$ to zero, i.e., $`W_{j,k}=0,`$ for $`j=1,2,\mathrm{},j_m`$.

Step 4:

In the fourth step, we reconstruct the denoised data $`\widehat{x}_d(t)`$ by taking inverse transform of the coefficients $`W_{j,k}`$ :

$`\widehat{x}_d(t)=c_\psi {\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}W_{j,k}\psi _{j,k}(t),`$ (4)
where $`c_\psi `$ is normalization constant given by

$`c_\psi =\mathrm{\hspace{0.33em}\hspace{0.33em}1}/_{\mathrm{}}^{\mathrm{}}\frac{|\widehat{\psi }(\omega )|^2}{\omega }𝑑\omega <\mathrm{},`$

with $`\widehat{\psi }(\omega )`$ as the Fourier transform of the wavelet function $`\psi (t)`$.

Step 5:

In the fifth and final step $`\widehat{x}_d(t)`$ is integrated to yield the cleansed signal $`\widehat{x}(t)`$:

$`\widehat{x}(t)={\displaystyle \widehat{x}_d(t)𝑑t}.`$ (5)
There exists a commutativity property between the operation of differentiation/integration and wavelet transform. Therefore first differentiating the signal and then taking DWT is equivalent to carrying out the two operations in reverse order. This implies that the same result can be obtained by switching the order between the first and second steps, and then between the fourth and fifth.

The effectiveness of the method lies in the following observations. Upon differentiation, contribution due to white noise moves towards the finer scales because the process of differentiation converts the uncorrelated stochastic process to a first order moving average process and thereby distributes more energy to the finer scales. That the differentiation of white noise brings about this behavior is known in the Fourier spectrum (Box et. al., 1994). It may be noted that the nature and effectiveness of separation depend on the wavelet basis function chosen and also on the properties of the derivatives of WT, which is in itself a highly interesting and not fully understood subject (Strang and Nguyen, 1996). For the signals studied in this paper, model as well as experimental, the wanted signal features lie in the coarser wavelet scales while the unwanted signal features after differentiation lie in the finer resolution wavelet scales. This is because the size of the data set handled decides the total number of scales available and a suitable choice can bring out the noisy signal WT features lying in the coarser scales. This justifies the assumption that fine scale features can be removed by setting the corresponding wavelet coefficients to zero and coarse scale features retained after differentiation. For this reason we also see a clear separation in the scales attributed to noise and those for the signal. A threshold scale for noise removal is thus identified and this leads to an automation for noise removal.

Result and Discussion

We first take up the three model systems and discuss the observations. For test purposes the pure signal obtained from these systems were corrupted with noise of certain strength. For the systems chosen for study, viz., Lorenz, Autocatalator, and Rössler, Table 1 summarizes the details, i.e. the equations governing their dynamics, the values chosen for the parameters, and the nature of the evolution of these systems for these sets of parameter values. These values were chosen appropriately so that the dynamics is chaotic. In our initial studies, purely for testing purposes, we studied situations where we ensure that all scales are affected by noise. In the wavelet domain this can be conveniently carried out by perturbing the wavelet coefficients in the following way. The differential equations are first numerically integrated to obtain pure signal $`x^0(t_i)`$ at equidistant time steps $`t_i`$. We then take DWT of the signal, and add white noise $`\eta `$ of zero mean and certain strength, i.e., $`W_{j,k}=W_{j,k}^0+\eta `$, where $`W_{j,k}^0,W_{j,k}`$ are the wavelet coefficients of pure and noisy signals respectively. We take the strength of noise as the relative percentage of the difference between the maximum and minimum of the signal value. Since each coefficient $`W_{j,k}`$ is individually affected by the noise, this procedure ensures equal weightage for presence of noise at all scales. Reconstructing the time series signal with this perturbed set of wavelet coefficients gave us the noisy signal to be cleansed. Our studies in this fashion did show that noise and signal separation was achieved. For the subsequent studies we followed the usual way of corrupting the signal by additive noise, i.e.,

$`x(t_i)=x^0(t_i)+\eta (t_i),`$ (6)
where $`\eta (t_i)[.5,.5]`$ is the noise with zero mean and uniform distribution, and $`J`$ the number of available dyadic scales. We have taken data size of 16384 ($`=2^{14}`$) points for all these three systems, and so $`J=14`$.

In Fig. 1 we plot our observations for the Lorenz system. Fig. 1 (a) shows the power at different scales in the pure signal $`x^0(t_i)`$. In Fig. 1 (b) we plot the scalewise power distribution after numerically differentiating the pure signal. We see that almost the entire power of the differentiated data is accumulated within the dyadic scales 4 and 9 (the signal power between scales 10 and 14 has disappeared by the process of differentiation). Fig. 1 (c) plots the scalewise power in the noisy signal $`x(t_i)`$ when the pure signal is infected with noise $`\eta (t_i)`$ of a typical strength of 5% of the signal (that is, 5% of the difference in the maximum and minimum values of $`x^0(t_i)`$). Because of the relative larger contribution at all scales from the pure signal, compared to that from noise, it is impossible to distinguish between the two components, and the figure looks qualitatively very similar to Fig. 1 (a). However, when we plot the scalewise power distribution of the differentiated noisy data in Fig. 1 (d), the signal contribution can easily be identified and also compared with the plot in Fig. 1 (b). It is to be noted that the difference in the values of the two peaks in Figs. 1 (b) and (d) arises because of the power being normalized by the respective total signal power. The contribution due to noise shows up in the finer scales. A clear minimum with close to zero value separates out two distinct regions. Fig. 1 (e) exhibits a small segment of the signal after the noise has been successfully removed following the procedure outlined above. All the three signals – pure, noisy, and cleansed – are overlaid for the sake of comparison.

In Fig. 2 we show the results for the Autocatalator and the Rössler reacting systems. Fig. 2 (a) and (c) plot, respectively, the scalewise power distribution for noise–infected signals obtained from the two systems. Like in the case of Lorenz system, it is evident that here also one cannot distinguish the noise and signal components. Fig. 2 (b) and (d) exhibit scalewise power profile for the differentiated data of the two signals respectively. The clear separation is again obvious.

In order to quantitatively estimate the efficiency of our denoising method, we have made the following error estimation (Kostelich and Schreiber, 1993) for the above three model systems. Since in all these cases the pure signal is known, a measure of the amount of error present in the cleaned data is obtained by taking rms deviation of the cleaned signal $`\widehat{x}(t_i)`$ from the pure signal $`x^0(t_i)`$ as follows,

$`\widehat{E}=\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\widehat{x}(t_i)x^0(t_i))^2\right)^{1/2},`$ (7)
where $`N`$ is the length of the time series. Similar quantity $`E`$ for the noisy data $`x(t_i)`$ is also computed. The condition $`\widehat{E}/E<1`$ guarantees that noise has been successfully reduced. The error estimator $`\widehat{E}/E`$ is a natural measure for noise reduction when the original dynamics is known (Kostelich and Schreiber, 1993). In Fig. 3 we plot $`\widehat{E}/E`$ against noise strength, for the three model systems. We see that for the entire range of noise values, and even with noise level as high as 10% of the signal exhibiting chaotic dynamics, $`\widehat{E}/E`$ remains appreciably below unity. Thus the plot demonstrates the efficiency of the approach. Different wavelet basis functions may change the nature and also improve the efficiency further.

We now discuss our method when applied to raw data obtained from two real chemical systems. In the first system, the time series data was obtained from the measurements of the pressure fluctuations in a fluidized bed, which consists of a vertical chamber inside of which a bed of solid particles is supported by an upwardly moving gas. Our system used a bed of silica sand particles (of mean diameter 200 microns) with a settled height of 500 mm, fluidized by ambient air in a transparent vessel 430 mm across and 15 mm wide. Beyond a critical inlet gas velocity, viz. the minimum bubbling velocity, the gas passes through the bed in the form of bubbles, thereby churning the solid and gas mixture in a turbulent manner. The time series data have been taken by measuring the pressure fluctuations inside this mixture, relative to atmospheric pressure, using a pressure transducer attached to a probe inserted into the fluid bed. The bed was operated at an inlet gas velocity of 0.85 m/sec, and the pressure fluctuations were recorded at a sampling rate of 333 Hz (333 data points per second). As a standard procedure, we normalize the data by subtracting mean and dividing by standard deviation (Constantinides, 1987; Bai et. al., 1997). In Fig. 4 we show the results obtained after the data have been subjected to denoising. Fig. 4 (a) shows the power distribution at different scales in the original experimental signal, while in Fig. 4 (b) we plot the scalewise power profile of the differentiated data. Again one clearly sees the two distinct contributions due to the noise and signal components. Fig. 4 (c) shows short segments of the denoised signal and the original signal which is overlaid for comparison. The cleaned signal is seen to be smooth indicating that the noise has been removed.

In the second chemical system, the time series data was obtained by sampling a measure related to the conductivity in a 3 liter liquid surfactant manufacturing experiment, at a sampling rate of 500 Hz. The time series is highly nonstationary since at various stages the operational parameters are altered (increasing the temperature for certain duration, then adding actives to the liquid, etc.). We studied unfiltered noisy data sets from the experiments, in order to check if our method can filter the noise out and also bring forth some intrinsic features of the system. We used our denoising algorithm to treat this data set in a slightly different way. The aim was to remove the finer scales from the differentiated data one by one, starting from the lowest (dyadic scale 1) and gradually going up, so that at each stage (after integrating the data) the observable frequencies in the filtered signal may be related to identifiable physical sources. Fig. 5 (a) shows a small segment (1 second long) of the noisy data. In Fig. 5 (b) we plot, on the same scale as in the earlier figure, the filtered data, using our method to remove the lowest dyadic scale 1. One can now clearly identify a 50 Hz component, due to the signal from electrical power supply (the ‘net frequency’). By removing scale 2 alongwith scale 1, the net frequency goes away, and the filtered data exhibits a 13 Hz component superimposed with occasional spikes. This 13 Hz signal shows up clearly in the filtered data with scale 3 also removed. The same Fig. 5 (b) shows this data, overlaid on the data with scale 1 removed. This 13 Hz may have arisen from the stirring device which has two blades and revolves with 260 rpm, coresponding to approximately 10 Hz. The electronic signal had an antialiasing feature of no more than 250 Hz and therefore aliasing (beating) may be ruled out. It may also be mentioned here that the Fourier power spectrum of the denoised signal shows a spike at 50 Hz frequency, whose removal resulted in a residual spectrum consisting mainly of a background continuum without any appreciable peak around 13 Hz. This study with the present example suggests that the wavelet transform methodology offers considerable benefits in the recovery of intrinsic signal components.

Summary

We have presented a new and alternative algorithm for noise reduction using discrete wavelet transform. We believe that our algorithm will be beneficial in various noise reduction applications, and that it shows promise in developing techniques which can resolve an observed signal into its various intrinsic components. In our method the threshold for reducing noise comes out automatically. The algorithm has been applied to three model flow systems - Lorenz, Autocatalator, and Rössler systems - all evolving chaotically. The method is seen to work quite well for a wide range of noise strengths, even as large as 10% of the signal level. We have also applied the method successfully to noisy time series data obtained from the measurement of pressure fluctuations in a fluidized bed, and also to that obtained by conductivity measurement in a liquid surfactant experiment. In all the illustrations we have been able to observe that there is a clean separation in the frequencies covered by the differentiated signal and white noise. However, if the noise is colored, a certain degree of overlap between the signal and noise may exist even after differentiation. For this complex situation, the method needs to be improved upon.

Acknowledgement

Authors acknowledge Unilever Research, Port Sunlight, for financial and other assistance. Part of the work has been carried out under the aegis of Indo–Australian S&T program DST/INT/AUS/I–94/97.

Literature cited

Abarbanel, H. D. I., “The Observance of Chaotic Data in Physical Syatems”, Rev. Mod. Phys. 65, 1340 (1993)

Bai, D., T. Bi, and J. R. Grace “Chaotic Behavior of Fluidized Beds Based on Pressure and Voidage Fluctuations”, AIChE Journal 43, 1357 (1997)

Bakshi, B., and G. Stephanopoulos “Representation of Process Trends. Part IV. Induction of Real–time Patterns from Operating Data for Diagnoses and Supervisory Control”, Computers chem. Engng. 18, 303 (1994)

Box, G. E. P., G. M. Jenkins and G. C. Reinsel, “Time Series Analysis Forecasting and Control” Prentice Hall, Englewood Cliffs (1994)

Carrier, J. F., and G. Stephanopoulos “Wavelet–Based Modulation in Control–Relevant Process Identification”, AIChE Journal 44, 341 (1998)

Cawley, R., and G.–H. Hsu “Local–Geometric–Projection Method for Noise Reduction in Chaotic Maps and Flows”, Phys. Rev. A 46, 3057 (1992)

Cohen, L., “Time frequency analysis” Prentice Hall, Englewood Cliffs (1995)

Constantinides, A., “Applied Numerical Methods with Personal Computers” McGraw–Hill Book Company, USA (1987)

Daubechies, I., “Ten Lectures on Wavelets” SIAM, Philadelphia (1990)

Donoho, D. L., and I. M. Johnstone “Adapting to Unknown Smoothness via Wavelet Shrinkage”, J. American Statistical Association 90, 1200 (1995)

Farmer, J. D., and J. J. Sidorowich “Optimal Shadowing and Noise Reduction”, Physica D 47, 373 (1991)

Grassberger, P., T. Schreiber, and C. Schaffrath “Non–linear Time Sequence Analysis”, Int. J. Bifurcation Chaos 1, 521 (1991)

Guckenheimer, E. and P. Holmes, “Nonlinear Oscillation, Dynamical Systems and Bifurcations of Vector Fields” Springer Verlag, Berlin (1983)

Härdle, W., “Applied Nonparametric Regression” Econometric Society Monographs, Cambridge University Press (1990)

Holschneider, M., “Wavelets: An Analysis Tool” Clarendon Press, Oxford (1995)

Kantz, H. and T. Schreiber, “Nonlinear Time Series Analysis” Cambridge University Press, Cambridge (1997)

Kostelich, E. J., and J. A. Yorke “Noise Reduction in Dynamical Systems”, Phys. Rev. A 38, 1649 (1988)

Lorenz, E. “Deterministic Nonperiodic Flow”, J. Atmos. Sc. 20, 130 (1963)

Luo, R., M. Misra, S. J. Qin, R. Barton, and D. M. Himmelblau “Sensor Fault Detection via Multiscale Analysis and Nonparametric Statistical Inference”, Ind. Eng. Chem. Res. 37, 1024 (1998)

Lynch, D. T. “Chaotic Behavior of Reaction Systems: Mixed Cubic and Quadratic Autocatalators”, Chem. Engg. Sci. 47, 4435 (1992)

Nason, G. P., “Wavelet Regression by Cross–validation” Dept. of Mathematics, University of Bristol (1994)

Press, W. H., B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, “Numerical Recipes” Cambridge University Press, Cambridge (1987)

Rössler, O. E. “Chaotic Behavior in Simple Reaction Syatems”, Z. Naturforsch. 31 a, 259 (1976)

Safavi, A. A., J. Chen, and J. A. Romagnoli “Wavelet–Based Density Estimation and Application to Process Monitoring”, AIChE Journal 43, 1227 (1997)

Sauer, T. “A Noise Reduction Method for Signals from Nonlinear Systems”, Physica D 58, 193 (1992)

Strang, G. and T. Nguyen, “Wavelets and Filter Banks” Wellesley–Cambridge Press (1996)

Strogatz, S. H., “Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering” Addison–Wesley (1994)

Figure Captions

* Plots for Lorenz system, with parameter values as stated in Table 1. 16384 ($`=2^{14}`$) data points are considered. Scalewise power distribution is plotted against the dyadic scale, (a) for the pure signal (the interpolated line through the data points is drawn for visualization), (b) using the data obtained after differentiating the signal, (c) for the signal corrupted by noise, and (d) using differentiated data of noisy signal. A segment of cleansed signal is shown in (e), alongwith the pure and noisy signals overlaid for comparison.
* Plots for Autocatalator and Rössler systems, in (a), (b) and (c), (d) respectively, for parameter values as in Table 1. Scalewise power profile plotted, (a) for the noisy autocatalytic signal, (b) using data after differentiating the signal, (c) for noisy Rössler signal, and (d) using the differentiated data of noisy signal.
* The error estimator $`\widehat{E}/E`$ plotted against the noise strength for all the three systems.
* Plots for the fluidized bed experiment. (a) Scalewise power profile is shown, for (a) the experimental signal, and (b) the data after the signal has been numerically differentiated. A small segment of the cleansed signal is shown in (c), alongwith the original signal for comparison.
* Plots for the liquid surfactant experiment. (a) A segment of the original noisy data, 1 second long. (b) The filtered data, with scale 1 removed, and with scales 1, 2 and 3 removed (on the same axes–scales as (a)).