text:
# Order statistics of Rosenstock’s trapping problem in disordered media

## I Introduction

Rosenstock’s trapping problem is a fundamental problem of random walk theory that has been of interest for decades TrappingN1 ; BluKlafZum ; Hughes ; Havlin . Most studies refer to the case in which a single ($`N=1`$) random walker is placed initially at a site of a Euclidean or disordered lattice which is randomly filled with trap sites at a concentration $`c`$, and then performs a random walk until it is absorbed by one of these traps. The statistical quantity of main interest in this problem is the survival probability $`\mathrm{\Phi }_1(t)`$ that the random walker is not trapped by time $`t`$, from which one can obtain the moments $`t^m=_0^{\mathrm{}}t^m[1\mathrm{\Phi }_1(t)]𝑑t`$ of the lifetime of this random walker. This problem has its origin in Smoluchowski’s theory of coagulation of colloidal particles Hughes ; Havlin ; Hollander  and has been applied to many systems in physics and chemistry such as trapping of mobile defects in crystals with point sinks Beeler ; Rosenstock ; Damask , the kinetics of luminescent organic materials Rosenstock , anchoring of polymers by chemically active sites Oshanin , and atomic diffusion in glasslike materials Miyagawa , among others.

A generalization of the trapping problem to the case of $`N`$ independent random walkers was studied by Krapivsky and Redner in Ref. KR . In particular, they studied the problem of $`N`$ diffusing predators placed initially at a given distance from a diffusing prey in one dimension. The model in which static preys are stochastically distributed all to one side of the predators was the subject of a later analysis OneSided . Shortly after, the order statistics of the trapping problem in $`d`$-dimensional Euclidean lattices for a set of $`N1`$ independent random walkers, i.e., the statistical description of the time $`t_{j,N}`$ elapsed until the first $`j`$ out of $`N`$ independent random walkers (initially starting at the same site) are trapped by quenched traps randomly arranged on $`d`$-dimensional Euclidean lattices, was studied (and rigorously solved for the one-dimensional case) in Ref. PREost . In this work, the moments $`t_{j,N}^m`$, $`m=1,2,3,\mathrm{}`$ were calculated from the probability $`\mathrm{\Phi }_{j,N}(t)`$ that $`j`$ random walkers of the initial set of $`N`$ have been absorbed by time $`t`$. The key step in this calculation was the assumption of independency of the random walkers that allowed the establishment of a relationship between $`\mathrm{\Phi }_{j,N}(t)`$ and the survival probability $`\mathrm{\Phi }_N(t)\mathrm{\Phi }_{0,N}(t)`$ of the full set of $`N`$ random walkers PREost . The survival probability $`\mathrm{\Phi }_N(t)`$ was calculated by means of Rosenstock’s approximation which required the evaluation of the first moment of the number $`S_N(t)`$ of different sites visited (territory explored) by $`N`$ random walkers.

Interest in multiparticle diffusion problems has had a boost lately because of some advances in optical spectroscopy SingMol  that make it possible to monitor events corresponding to single particles of an ensemble. The simultaneous tracking of $`N1`$ fluorescently labelled particles and the analysis of the motions of the particles allows the study of local conditions (mechanical response, visco-elasticity) inside many complex structures such as fibrous polymers and the intracellular medium MultiTracking . But biological samples and many real inorganic substances are disordered media (as opposed to translationally invariant Euclidean media) which are usually described as stochastic fractals Havlin ; FractalBooks ; Pfeifer . There are two main reasons for this identification: disordered systems share the statistical fractal structure of stochastic fractal models and diffusion is anomalous in both media.

The single-particle ($`N=1`$) Rosenstock trapping problem in fractal media have been thoroughly discussed by Blumen, Klafter and Zumofen BluKlafZum ; BluKlafZumPRB . In this paper we will study its *multiparticle*  ($`N1`$) version, which is relevant for all those cases where the diffusing particles are placed (or created) in bunches. This may be especially important if the first or first few particles that are absorbed lead to a trigger effect. Here, we generalize to disordered fractal substrates the results for the order statistics of the multiparticle trapping problem obtained for Euclidean media in Ref. PREost . We will discover that for the two-dimensional incipient percolation aggregate, and in sharp contrast with the Euclidean media results, the zeroth-order Rosenstock approximation is quite incapable of describing the survival probabilities, and therefore the order statistics of the multiparticle trapping problem we are dealing with. This is because the ratio between the cumulants $`\kappa _m`$ of the distribution of $`S_N(t)`$ and $`S_N(t)^m`$ is relatively large (and almost constant) for $`N1`$. We traced the origin of this behavior to the fact that the fluctuations in the number $`S_N(t)`$ of distinct sites explored by a large number $`N`$ of random walkers are negligible relative to the fluctuations in the number of sites that form the stochastic substrate. As a practical consequence, higher-order Rosenstock approximations are required for these media in order to give an account of the order statistics trapping problem with an accuracy similar to those reached by means of the zeroth-order Rosenstock approximation for Euclidean lattices. The idea of evaluating the survival probability for the multiparticle trapping problem for Euclidean media by means of the Rosenstock approximation was first suggested by Larralde *et al.*  in LarraldePRA , although, to the best of our knowledge, it has not been implemented (except for the so-called “one-sided trapping problem” OneSided ) perhaps for the lack of precise expressions for the moments $`S_N^m(t)`$ of $`S_N(t)`$. However, for the percolation aggregate, we will discover in Sec. III that one can get a very good estimate of $`S_N^m(t)`$ from the value of the first moment $`S_N(t)`$. This is fortunate because the asymptotic expansion of $`S_N(t)`$ for large $`N`$ is known for this medium SNtFrac .

The multiparticle Rosenstock trapping problem we are considering here can be seen as a stochastic generalization of the problem of the order statistics for the sequence of trapping times (or exit times) of a set of $`N`$ independent random walkers, all starting form the same site at the same time, when the traps form a “spherical” absorbing boundary with a fixed radius. This problem was first studied by Lindenberg et al.  Lindenberg  and Weiss et al.  WeisShulLind  for Euclidean lattices (mainly for the one dimensional case). After these pioneering works, improved results and extensions to deterministic and random fractal substrates have been reported trapesfera ; DragerKlafter . A closely related multiparticle Brownian problem has been recently considered by Bénichou *et al.* BenichouJPA03 : they studied the join residence times of $`N`$ independent Brownian particles in a disc of a given radius. In particular, they studied the time spent by all $`N`$ particles simultaneously in the disc within a given time interval, and the time which at least $`Nj`$ out of $`N`$ particles spend together in the disc within a time interval.

The plan of the paper is as follows. In Sec. II the expressions that describe the order statistics of the trapping process are deduced. In Sec. III we study the moments $`S_N^m(t)`$ of the territory explored by $`N`$ independent random walkers on two-dimensional incipient percolation aggregates by means of numerical simulation. The results of Sec. III are applied in Sec. IV to obtain the survival probability $`\mathrm{\Phi }_N(t)`$ by means of Rosenstock’s approximation. Then we calculate the moments $`t_{j,N}^m`$ of the time elapsed until the first $`j`$ random walkers are trapped for every $`j=1,2,\mathrm{}`$ and $`m=1,2,\mathrm{}`$ and compare these predictions with simulation results for the two-dimensional incipient percolation aggregate. A general discussion and conclusions are given in Sec. V.

## II Definitions and fundamental relations

The results and definitions of this section have already been discussed in detail in the context of the trapping problem in Euclidean media PREost . However, we will briefly summarize those results that are basic and necessary in order to follow the arguments in the rest of the paper.

Let us first define $`\mathrm{\Psi }_{j,N}(t)`$ as the probability that $`j`$ random walkers of the initial set of $`N`$ have been absorbed by time $`t`$ by a given configuration of traps arranged on a given realization of the disordered substrate. The quantity of statistical interest is the average $`\mathrm{\Phi }_{j,N}(t)=\mathrm{\Psi }_{j,N}(t)`$ performed over all the possible outcomes of the “trapping experiment” carried out in a quenched configuration of traps in a given lattice realization followed by an average over all trap configurations and lattice realizations. We will also denote by $`\mathrm{\Psi }(t)`$ the probability that a single random walker has not been absorbed by time $`t`$ in this quenched configuration of traps placed upon a specific lattice realization. This is commonly known as the survival probability. It is then clear that

$$\mathrm{\Psi }_{j,N}(t)=\left(\genfrac{}{}{0pt}{}{N}{j}\right)\left(1\mathrm{\Psi }\right)^j\mathrm{\Psi }^{Nj}=\left(\genfrac{}{}{0pt}{}{N}{j}\right)\underset{m=0}{\overset{j}{}}(1)^m\left(\genfrac{}{}{0pt}{}{j}{m}\right)\mathrm{\Psi }^{Nj+m}=(1)^j\left(\genfrac{}{}{0pt}{}{N}{j}\right)^j\mathrm{\Psi }_{0,N}$$
(1)
where $`^j\mathrm{\Psi }_{0,N}(t)=_{m=0}^j(1)^m\left(\genfrac{}{}{0pt}{}{j}{m}\right)\mathrm{\Psi }_{0,Nm}(t)`$ is just the backward difference formula for the $`j`$th derivative of $`\mathrm{\Psi }_{0,N}(t)`$, $`d^j\mathrm{\Psi }_{0,N}(t)/dN^j`$. Averaging over different configurations, and taking into account that $`\mathrm{\Phi }_{0,N}(t)\mathrm{\Phi }_N(t)=\mathrm{\Psi }_{0,N}(t)`$ and $`\mathrm{\Phi }_{j,N}(t)=\mathrm{\Psi }_{j,N}(t)`$, we get

$$\mathrm{\Phi }_{j,N}(t)=(1)^j\left(\genfrac{}{}{0pt}{}{N}{j}\right)^j\mathrm{\Phi }_N(t).$$
(2)
Let us call $`𝒮_j`$ the state in which $`j`$ particles have been absorbed and $`Nj`$ particles of the initial set of $`N`$ are still diffusing, and let $`h_{j,N}(t)dt`$ be the probability that the $`j`$th absorbed particle of the initial set of $`N`$ disappears during the time interval $`(t,t+dt]`$. The change of the probability of $`𝒮_j`$ during the time interval $`(t,t+dt]`$ is given by $`\mathrm{\Phi }_{j,N}(t+dt)\mathrm{\Phi }_{j,N}(t)`$. But this probability changes during this time interval by two causes: first, because the state $`𝒮_{j1}`$ can become the state $`𝒮_j`$ if a particle of the set of $`Nj+1`$ particles still diffusing is trapped during the time interval $`(t,t+dt]`$ \[which happens with probability $`h_{j,N}(t)dt`$\] and, second, because $`𝒮_j`$ can become the state $`𝒮_{j+1}`$ if a particle of the $`Nj`$ particles still diffusing is trapped during this time interval \[which happens with probability $`h_{j+1,N}(t)dt`$\]. Therefore $`\mathrm{\Phi }_{j,N}(t+dt)\mathrm{\Phi }_{j,N}(t)=[h_{j,N}(t)h_{j+1,N}(t)]dt`$, i.e.,

$$h_{j+1,N}(t)=h_{j,N}\frac{d}{dt}\mathrm{\Phi }_{j,N}(t)$$
(3)
with $`h_{0,N}=0`$. Then, the $`m`$th moment of the time in which the $`j`$th particle is trapped is given by

$$t_{j,N}^m=_0^{\mathrm{}}t^mh_{j,N}(t)𝑑t.$$
(4)
Using Eqs. (2), (3) and (4) we find the *exact* recursion relation

$$t_{j+1,N}^m=t_{j,N}^m+(1)^j\left(\genfrac{}{}{0pt}{}{N}{j}\right)^jt_{1,N}^m$$
(5)
where

$$t_{1,N}^m=m_0^{\mathrm{}}t^{m1}\mathrm{\Phi }_N(t)𝑑t.$$
(6)
The set of Eqs. (5) and (6) is remarkable because it implies that the order statistic of the trapping problem can be described from the knowledge of $`t_{1,N}^m`$ only. The difference derivative $`^j`$ in Eqs. (2) and (5) can be approximated by the ordinary derivative $`d^j/dN^j`$ when $`jN`$. This will be justified (and used) in Sec. IV.

## III Moments of the territory explored by $`N1`$ random walkers on a two-dimensional incipient percolation aggregate

The diffusion in percolation clusters as a model of transport in disordered media was first suggested by de Gennes Gennes . Percolation clusters are disordered fractals: they share the self-similarity property with deterministic fractals build up through deterministic rules but only in a statistical sense. In order to characterize these fractals several static and dynamic exponents have been defined. Perhaps, the most widely known is the fractal dimension, $`d_f`$, which, in the case of disordered systems, is more conveniently defined using the scaling of mass with linear size, $`ML^{d_f}`$. However, random walkers in disordered structures are forced to follow the paths formed by the bonds between sites and, consequently, it is more natural to define a chemical (or topological) distance between two sites as the length of the shortest path along lattice bonds, $`\mathrm{}`$. If we consider the number of sites inside an hypersphere of radius $`\mathrm{}`$, $`V(\mathrm{})`$, usually known as chemical volume (also coincides with the mass if we assume that every site has a unit mass) it is expected that $`V(\mathrm{})V_0\mathrm{}^d_{\mathrm{}}`$, where the brackets refers to an average over all possible realization of the lattices and $`d_{\mathrm{}}`$ is the chemical dimension. Similarly, the generalized Einstein’s law of diffusion for anomalous systems can be written in terms of the ordinary Euclidean distance, $`r`$, or the chemical one, $`\mathrm{}`$, and we have:

$$\begin{array}{ccc}\hfill r^2& & 2Dt^{2/d_w},\hfill \\ \multicolumn{3}{c}{}\\ \hfill \mathrm{}^2& & 2D_{\mathrm{}}t^{2/d_w^{\mathrm{}}},\hfill \end{array}$$
(7)
for $`t1`$, and where $`D`$ and $`D_{\mathrm{}}`$ are the diffusion coefficient and the chemical diffusion coefficient, respectively. The exponent $`d_w`$ is the random walk dimension, also known as diffusion exponent. The exponent $`d_w^{\mathrm{}}`$ corresponding to the chemical (or topological) metric is called chemical random walk dimension. Another important exponent appearing in the theory of random walks in disordered media is the spectral or fracton dimension, $`d_s=2d_f/d_w=2d_{\mathrm{}}/d_w^{\mathrm{}}`$ Havlin .

Rosenstock’s procedure Rosenstock  for evaluating the survival probability of a set of random walkers requires the knowledge of the first moments of the territory explored $`S_N(t)`$ by these random walkers. This is an interesting (and difficult) problem in itself that has already been thoroughly studied in the case of $`N1`$ independent random walkers KR ; LarraldePRA ; LarraldeNPO ; DragerKlafter ; SNtEuc ; SNtFrac ; MultiparticleReview  although only the first moment $`S_N(t)`$ has been rigourously estimated SNtEuc ; SNtFrac ; MultiparticleReview . The average value of the territory explored by $`N1`$ random walkers, all starting from the same site, in a disordered medium was analyzed in Ref. DragerKlafter ; SNtFrac  and it was found that SNtFrac

$$S_N(t)\overline{S}_Nt^{d_s/2}$$
(8)
with

$$\overline{S}_N=V_0(2D_{\mathrm{}})^{d_{\mathrm{}}/2}\left(\frac{\mathrm{ln}N}{\widehat{c}}\right)^{d_{\mathrm{}}/v}\left[1d_{\mathrm{}}\frac{d_w^{\mathrm{}}1}{d_w^{\mathrm{}}}\underset{n=1}{\overset{\mathrm{}}{}}\left(\mathrm{ln}N\right)^n\underset{m=0}{\overset{n}{}}s_m^{(n)}\left(\mathrm{ln}\mathrm{ln}N\right)^m\right].$$
(9)
The parameters $`\widehat{c}`$ and $`s_m^{(n)}`$ are characteristic of the lattice and some of their values for several Euclidean and fractal media are known SNtEuc ; SNtFrac ; MultiparticleReview . In particular, for the two-dimensional incipient percolation aggregate, Monte Carlo simulations in this substrate (with particles jumping from a site to one of its nearest neighbors placed at one unit distance in each unit time) have shown SNtFrac ; MultiparticleReview  that $`d_{\mathrm{}}1.65`$, $`v=d_w^{\mathrm{}}/(d_w^{\mathrm{}}1)`$ with $`d_w^{\mathrm{}}2.45`$, $`d_s=2d_{\mathrm{}}/d_w^{\mathrm{}}1.35`$, $`V_01.1`$, $`\widehat{c}1.05`$, $`2D_{\mathrm{}}1.2`$, $`s_0^{(1)}=\gamma \mathrm{ln}\widehat{A}\widehat{c}^{\widehat{\mu }}0.62`$ ($`\gamma 0.577216`$ is the Euler constant) and $`s_1^{(1)}=\widehat{\mu }0.8`$ (see Table 1). Hence we have a reasonable estimate of the asymptotic series for $`S_N(t)`$ in Eq. (9) up to first order ($`n=1`$), which is sufficient to account for simulation results, as Fig. 1 shows.

In our simulations we also evaluated the second cumulant (variance), $`\kappa _2(S_N)=S_N^2S_N^2`$, and the third cumulant, $`\kappa _3(S_N)=S_N^33S_N^2S_N+2S_N^3`$, of the territory explored $`S_N(t)`$ as they are necessary for implementing the extended Rosenstock approximation (see Sec. IV). We found that the ratio $`k_m\kappa _m(S_N)/S_N(t)^m`$, although not very sensitive to the value of $`N`$ (one can see in Figs. 2 and 3 that these parameters are well represented by $`k_2=0.13\pm 0.02`$ and $`k_3=0.015\pm 0.02`$ over a wide range of $`N`$ values) seems to tend to a constant value for large $`N`$ (about $`0.14`$ for $`k_2`$ and $`0.016`$ for $`k_3`$). This is a surprising behavior that departs considerably from that of Euclidean media. For example, for the $`d`$-dimensional Euclidean lattices it was found that $`k_2`$ goes as $`1/\mathrm{ln}^2N`$ for large $`N`$. Also, the value of $`k_2`$ for large $`N`$ is much smaller than for the percolation aggregate (for example, for $`N=2^{10}`$, $`k_21/15^2,1/30^2`$ and $`1/50^2`$ in the one-, two- and three-dimensional Euclidean lattices, respectively PREost ), which has important consequences for the accuracy of Rosenstock’s approximations of different orders, as we will show in Sec. IV. The disordered nature of the substrate must be the reason for these remarkable differences in the behavior of $`k_m`$. What is happening is that, for large $`N`$, the fluctuations in the number $`S_N(t)`$ of distinct sites explored by a large number $`N`$ of random walkers are dominated by the fluctuations (over the set of stochastic lattice realizations used in the simulations) of the number of sites inside a hypersphere of chemical radius $`\mathrm{}=\sqrt{2D_{\mathrm{}}}t^{1/d_w^{\mathrm{}}}`$. We summarize this claim in a conjecture as follows

$$\underset{N\mathrm{}}{lim}k_m=\underset{N\mathrm{}}{lim}\frac{\kappa _m(S_N(t))}{S_N(t)^m}=\frac{\kappa _m(V)}{V(\mathrm{})^m}\mathrm{}_m,$$
(10)
where $`V(\mathrm{})`$ is the chemical volume (number of sites) of a hypersphere of chemical radius $`\mathrm{}`$ and $`\kappa _m(V)`$ is the $`m`$th-cumulant of the distribution of $`V`$. Rigorously, the distance $`\mathrm{}`$ appearing in Eq. (10) is given by $`\sqrt{2D_{\mathrm{}}}t^{1/d_w^{\mathrm{}}}`$ which is the radius of the diffusion front in the thermodynamic limit ($`N\mathrm{}`$). However, the quotient $`\mathrm{}_m`$ is not very sensitive to $`\mathrm{}`$ if a sufficiently large value of $`\mathrm{}`$ is taken. In Fig. 4 we plot a histogram for the chemical volume of a two-dimensional incipient percolation aggregate with $`\mathrm{}=100`$, evaluated using $`2000`$ realizations of the lattice. Thereby we find that $`\mathrm{}_20.14`$ and $`\mathrm{}_30.015`$ in very good agreement with the values of $`k_2`$ and $`k_3`$, respectively, for large $`N`$ (see figures 2 and 3). Consequently, we conclude that the fluctuations in $`S_N(t)`$ are dominated by the disorder of the substrate and the influence of the value of $`N`$ is completely overshadowed. Similar arguments were presented by Rammal and Toulouse in their pioneer work on percolation clusters RT .

## IV Order statistics of the trapping process

Assume that one has a quenched configuration of traps randomly placed on a given realization of the disordered lattice with probability $`c`$. If $`N`$ random walkers start from an origin site free from traps at $`t=0`$, it is clear that the probability that all random walkers survive by time $`t`$ is given by $`(1c)^{S_N(t)}`$. The average of this quantity over all possible random walks, trap configurations, and substrate realizations is known as the survival probability: $`\mathrm{\Phi }_N(t)=(1c)^{S_N(t)}`$. Using a well-known theorem in statistics, we can define the $`n`$th-order Rosenstock approximation $`\mathrm{\Phi }_N^{(n)}(t)`$ for estimating $`\mathrm{\Phi }_N(t)`$ as

$$\mathrm{\Phi }_N^{(n)}(t)=\mathrm{exp}\left[\underset{j=1}{\overset{n+1}{}}\frac{(\lambda )^j}{j!}\kappa _j(S_N)\right]$$
(11)
where $`\lambda \mathrm{ln}(1c)`$ and $`\kappa _j(S_N)`$ is the $`j`$th-cumulant of the distribution of the territory explored. In the limit $`n\mathrm{}`$ we recover the exact result for $`\mathrm{\Phi }_N(t)`$. In the case of the single particle ($`N=1`$) trapping problem, Eq. (11) is known as the extended Rosenstock approximation or truncated cumulant expansion Hughes ; Havlin ; Hollander ; BluKlafZumPRB ; ZumofenBlumen ; Weiss . Its generalization to the $`N`$ particle case was used in a one-dimensional trapping problem in Ref. OneSided .

In the previous section we showed that Monte Carlo simulations strongly suggest that $`\kappa _n(S_N)k_nS_N(t)^n`$ for large $`N`$, where $`k_n`$ are constants. \[Notice that $`\kappa _1(S_N)=S_N(t)`$, so that $`k_1=1`$.\] Therefore, inserting this result into Eq. (11), the $`n`$th-order Rosenstock approximation becomes

$$\mathrm{\Phi }_N^{(n)}(t)\mathrm{exp}\left[\underset{j=1}{\overset{n+1}{}}\frac{(\lambda )^j}{j!}k_jS_N^j\right]$$
(12)
or equivalently, by using Eq. (8),

$$\mathrm{\Phi }_N^{(n)}(t)\mathrm{exp}\left[\underset{j=1}{\overset{n+1}{}}\frac{(\lambda )^j}{j!}k_j\overline{S}_N^jt^{jd_s/2}\right].$$
(13)
We can now evaluate an approximation for the moments of the first trapping time, $`t_{1,N}^m`$, by means of Eq. (6) assuming that the contribution of $`\mathrm{\Phi }_N(t)`$ to $`t_{1,N}^m`$ is negligible for those times for which $`\mathrm{\Phi }_N(t)`$ and $`\mathrm{\Phi }_N^{(n)}(t)`$ differ substantially. Therefore, the substitution of Eq. (13) into Eq. (6) yields

$$t_{1,N}^mt_{1,N}^m_nm_0^{\mathrm{}}t^{m1}\mathrm{exp}\left[\underset{j=1}{\overset{n+1}{}}\frac{(\lambda )^j}{j!}k_j\overline{S}_N^jt^{jd_s/2}\right]𝑑t.$$
(14)
Writing $`v=\lambda \overline{S}_Nt^{d_s/2}`$, the $`n`$th-order Rosenstock’s estimate $`t_{1,N}^m_n`$ for $`t_{1,N}^m`$ becomes:

$$t_{1,N}^m_n=\frac{2m}{d_s}\left(\lambda \overline{S}_N\right)^{2m/d_s}\tau _n(m)$$
(15)
where

$$\tau _n(m)=_0^{\mathrm{}}v^{2m/d_s1}\mathrm{exp}\left[\underset{j=1}{\overset{n+1}{}}\frac{(1)^j}{j!}k_jv^j\right]𝑑v.$$
(16)
Therefore, we find that the different $`n`$th-order Rosentock approximations $`t_{1,N}^m_n`$ differ from each other only by a numerical factor $`\tau _n(m)`$ (an integral) that depends only on the substrate through its spectral dimension $`d_s`$ and the set of parameters $`k_j`$, $`j=1,2,\mathrm{}`$ that come from the distribution of the chemical volume of this substrate \[c.f. Eq. (10)\]. The integral in Eq. (16) is trivial for $`n=0`$ and yields $`\tau _0(m)=\mathrm{\Gamma }(2m/d_s)`$. Using the values of Table 1 we get $`d_s=2d_{\mathrm{}}/d_w^{\mathrm{}}1.375`$ and the estimates $`\tau _0(1)0.89`$, $`\tau _0(2)1.95`$ and $`\tau _0(3)10.9`$ for the two-dimensional incipient percolation aggregate. The integral in Eq. (16) only converges for even values of $`n`$ so the next meaningful approximation corresponds to $`n=2`$. Taking the values $`k_2=0.13`$ and $`k_3=0.015`$ (which describe $`k_2`$ and $`k_3`$ well over the range of values of $`N`$ used in our simulations: see Figs. 2 and 3) and evaluating the integral in Eq. (16) numerically, we found the second-order prefactors $`\tau _2(m)`$: $`\tau _2(1)1.24`$, $`\tau _2(2)6.04`$ and $`\tau _2(3)95.6`$, which are systematically much larger than the zeroth-order ones $`\tau _0(m)`$, especially when the order $`m`$ of the moment is large. This means that Rosenstock’s approximations of order higher than zero must be necessary to provide reasonable predictions for $`\mathrm{\Phi }_N(t)`$ and $`t_{1,N}^m`$ in disordered media, especially when the moment $`m`$ is large. It should be noticed that the expression for the first trapping time in Eq. (15) includes two approximations of different nature: (a) one due to the fact that we are using a finite number $`n`$ of terms in the cumulant expansion \[which only affects the factor $`\tau _n(m)`$\]; and (b) the other due to the finite number of terms considered for estimating $`\overline{S}_N`$ by means of the asymptotic series (9). Consequently, it is convenient to classify these approximations by indexing them with a pair of integers $`(n,l)`$: the first index gives the order of the Rosenstock approximation that is used, and the second gives the number of terms considered in the evaluation of $`S_N(t)`$. In this way, the approximation $`(n,0)`$ corresponds to the replacement in Eq. (15) of $`\overline{S}_N`$ by the leading term of the series of Eq. (9), so that $`t_{1,N}^m_n=t_{1,N}^m_{n0}[1+O(1/\mathrm{ln}N)]`$ with

$`t_{1,N}^m_{n0}`$ $`=`$ $`{\displaystyle \frac{2m\tau _n(m)}{d_s\left[\lambda V_0(2D_{\mathrm{}})^{d_{\mathrm{}}/2}\right]^{2m/d_s}}}\left({\displaystyle \frac{\mathrm{ln}N}{\widehat{c}}}\right)^{m(1d_w^{\mathrm{}})}`$ (17)
$`=`$ $`{\displaystyle \frac{T_n(m)}{\lambda ^{2m/d_s}(\mathrm{ln}N)^{m(d_w^{\mathrm{}}1)}}},`$ (18)
and where we have absorbed all the dependence on the lattice characteristic parameters ($`V_0`$, $`D_{\mathrm{}}`$, $`d_s`$, …) into the coefficient $`T_n(m)`$. In the same way, if we take the two first terms in the asymptotic series of of Eq. (9), we find $`t_{1,N}^m_n=t_{1,N}^m_{n1}[1+O(1/\mathrm{ln}N)^2]`$ where the approximation $`(n,1)`$ is

$$t_{1,N}^m_{n1}=\frac{T_n(m)}{\lambda ^{2m/d_s}(\mathrm{ln}N)^{m(d_w^{\mathrm{}}1)}}\left[1d_{\mathrm{}}\frac{d_w^{\mathrm{}}1}{d_w^{\mathrm{}}}\frac{s_0^{(1)}+s_1^{(1)}\mathrm{ln}N}{\mathrm{ln}N}\right]^{2m/d_s}$$
(19)
or, for $`\mathrm{ln}N1`$,

$$t_{1,N}^m_{n1}=\frac{T_n(m)}{\lambda ^{2m/d_s}(\mathrm{ln}N)^{m(d_w^{\mathrm{}}1)}}\left[1+m(d_w^{\mathrm{}}1)\frac{s_0^{(1)}+s_1^{(1)}\mathrm{ln}N}{\mathrm{ln}N}\right].$$
(20)
In Fig. (5) we compare simulation results for the trapping time of the first particle with the theoretical predictions given by Eqs. (18) and (19) when the parameters of Table 1 are used. We see that the second-order Rosenstock approximation leads to much better results than the standard zeroth-order approximation.

The moments $`t_{j,N}^m`$, $`j=2,3,\mathrm{}`$ corresponding to the trapping of the $`j`$th particle absorbed by the traps can be easily estimated by means of Eq. (5). However, we can also obtain an explicit expression for $`t_{j,N}^m`$ if we approximate the difference operator $`^j`$ in (5) by the derivative $`d^j/dN^j`$. The error in this approximation can be estimated from the equation

$$\frac{^jf(N)}{(\mathrm{\Delta }N)^j}=\frac{d^jf(N)}{dN^j}+O\left(\frac{d^{j+1}f(N)}{dN^{j+1}}\mathrm{\Delta }N\right).$$
(21)
In our case $`f(N)`$ is $`t_{1,N}^m`$ and $`\mathrm{\Delta }N=1`$ so that

$$^jt_{1,N}^m=\frac{d^jt_{1,N}^m}{dN^j}+O\left(\frac{d^{j+1}t_{1,N}^m}{dN^{j+1}}\right).$$
(22)
As

$$\frac{d^j}{dN^j}(\mathrm{ln}N)^\mu =(1)^j\mu \frac{(j1)!}{N^j}(\mathrm{ln}N)^{\mu 1}+\frac{\mathrm{ln}\mathrm{ln}N}{N^j}O[(\mathrm{ln}N)^{\mu 2}]$$
(23)
one finds, from Eq. (18) \[or Eq.(19)\], that

$$^jt_{1,N}^m=\frac{d^jt_{1,N}^m}{dN^j}\left[1+O(N^1)\right].$$
(24)
Taking into account that $`\left(\genfrac{}{}{0pt}{}{N}{j}\right)/N^j1/j!`$ for $`jN`$, we obtain from Eqs. (5), (18) \[or (19)\] and (23) the recursion relation

$$t_{j+1,N}^m=t_{j,N}^m+\frac{1}{j}m(d_w^{\mathrm{}}1)T_n(m)\lambda ^{2m/d_s}(\mathrm{ln}N)^{m(d_w^{\mathrm{}}1)1}\left[1+O\left(\frac{1}{\mathrm{ln}N}\right)\right]$$
(25)
which can be easily solved:

$$t_{j,N}^m=t_{1,N}^m+m(d_w^{\mathrm{}}1)T_n(m)\lambda ^{2m/d_s}\frac{\psi (j)+\gamma }{(\mathrm{ln}N)^{m(d_w^{\mathrm{}}1)+1}}\left[1+O\left(\frac{1}{\mathrm{ln}N}\right)\right]$$
(26)
where

$$\psi (j)=\psi (1)+\underset{r=1}{\overset{j1}{}}\frac{1}{r}$$
(27)
is the psi (digamma) function abramo , $`\psi (1)=\gamma `$, and $`\gamma `$ is the Euler constant. Equation (26) yields

$$t_{j,N}^m_{n0}=\frac{T_n(m)}{\lambda ^{2m/d_s}(\mathrm{ln}N)^{m(d_w^{\mathrm{}}1)}}$$
(28)
by using Eq. (18), and

$$t_{j,N}^m_{n1}=\frac{T_n(m)}{\lambda ^{2m/d_s}(\mathrm{ln}N)^{m(d_w^{\mathrm{}}1)}}\left[1+m(d_w^{\mathrm{}}1)\frac{\psi (j)+\gamma +s_0^{(1)}+s_1^{(1)}\mathrm{ln}N}{\mathrm{ln}N}\right]$$
(29)
when Eq. (20) is used. In Fig. (6) we compare the predictions for $`t_{2,N}`$ obtained from Eqs. (28) and (29) with simulation results. The results are similar to that found in Fig. (5) for $`t_{1,N}`$. In Fig. (7) the differences $`t_{j+1,N}t_{j,N}`$ estimated from Eq. (25) are also plotted in a scaled form for $`j=1,2`$. The theoretical prediction is that, for large $`N`$, these points should tend to lie along a straight line (which is true) with a slope $`(d_w^{\mathrm{}}1)T_2(1)/j`$, i.e., a slope $`0.76`$ for $`j=1`$ and 1.0 for $`j=2`$. The last prediction is not good for $`j=2`$ (the fitted value is 0.90), but this should not be surprising because in Eq. (25) we have ignored correction terms of order $`1/\mathrm{ln}N`$, which are very large even for huge values of $`N`$. The only way to remedy this deficiency would be by increasing the number of asymptotic terms retained in the evaluation of $`S_N(t)`$, which in turns requires knowing $`s_m^{(n)}`$ for $`n2`$ \[c.f. Eq. (9)\]. Unfortunately, these values of $`s_m^{(n)}`$ are very difficult to estimate by means of numerical simulations SNtFrac  and are unknown for $`n2`$.

The variance of $`t_{j,N}`$ is easily obtained from Eq. (15):

$$\sigma _{j,N}^2=\frac{[T_n(2)T_n^2(1)]}{\lambda ^{4/d_s}(\mathrm{ln}N)^{2(d_w^{\mathrm{}}1)}}\left[1+2(d_w^{\mathrm{}}1)\frac{\psi (j)+\gamma +s_0^{(1)}+s_1^{(1)}\mathrm{ln}N}{\mathrm{ln}N}+O\left(\frac{1}{\mathrm{ln}^2N}\right)\right]$$
(30)
and, consequently,

$$\frac{\sigma _{j,N}^2}{t_{1,N}^2}=d_s\frac{\tau _n(2)}{\tau _n(1)^2}1.$$
(31)
This is an interesting result because it means that the ratio between the variance of the first trapping time and the mean of that time is, for large $`N`$, independent of $`N`$. The numerical value of the ratio $`\sigma _{j,N}/t_{1,N}`$ is $`[d_s\mathrm{\Gamma }(4/d_s)/\mathrm{\Gamma }^2(2/d_s)]11.53`$ for the zeroth-order Rosenstock approximation, and $`2.0`$ for the second-order approximation. In Fig. 8 we plot this ratio for $`j=1,2,3`$ versus $`\mathrm{ln}N`$, and the second-order theoretical limit $`\sigma _{j,N}/t_{1,N}2`$ seems to be consistent with the simulation data.

Some considerations about the range of validity of the approximations developed in this paper are called for at this point. The approximation for $`S_N(t)`$ in Eqs. (8) and (9) is only valid in the so called Regime II or intermediate time regime SNtFrac ; MultiparticleReview . As the integral in Eq. (6) was evaluated assuming that the expression for $`S_N(t)`$ was valid for all times we conclude that the integral of $`mt^{m1}\mathrm{\Phi }_N(t)`$ over the short-time interval $`[0,t_\times ]`$ ($`t_\times \mathrm{ln}N`$ being the crossover time between Regime I and Regime II) has to be negligible relative to $`t_{1,N}^m`$, or equivalently $`(\mathrm{ln}N)^mt_{1,N}^m`$, for our approach and our results being reasonable. Taking into account the estimate for $`t_{1,N}^m`$ given in Eq. (18), this condition can be written as $`\lambda (\mathrm{ln}N)^d_{\mathrm{}}`$. The concentrations of traps we have used in our simulations, $`c=0.001`$ and $`c=0.008`$, verify this condition for all the values of $`N`$ considered. Apart from this upper bound on $`c`$, we must also point out that, as also for Euclidean lattices, our results break down if most of the trapping takes place within the long-time Donsker-Varadhan regime DV . Further reference to this limitation of the theory presented in this paper will be made below.

## V Summary and conclusions

We have dealt with the following order statistics problem: when $`N`$ independent random walkers all starting from the same site diffuse on a disordered lattice populated with a concentration $`c`$ of static trapping sites, what is the distribution of the elapsed times, $`t_{j,N}`$, until the first $`j`$ random walkers are trapped? We were able to generalize the theory developed for the special case of Euclidean lattices PREost  to the case of disordered substrates, and asymptotic expressions for the moments $`t_{j,N}^m`$ with $`jN`$ were obtained. To this end, we used the so-called Rosenstock approximation, which is suitable for not very large times and small concentrations of traps, $`c`$. In this approximation the survival probability of the full set of $`N`$ random walkers, $`\mathrm{\Phi }_N(t)`$, is expressed in terms of the cumulants $`\kappa (S_N)`$ of the distribution of the territory covered $`S_N(t)`$.

Monte Carlo simulation results for $`\kappa (S_N)`$ in the two-dimensional incipient percolation aggregate showed that for large $`N`$ the ratio $`k_m=\kappa _m(S_N)/S_N(t)^m`$ with $`m=2,3`$ hardly depends on $`N`$ and is very large in comparison with the corresponding Euclidean ratio. We attribute this behavior to the fluctuations in $`S_N(t)`$ being dominated by the fluctuations in the volume of the medium inside a hypersphere of chemical radius $`\mathrm{}\sqrt{2D_{\mathrm{}}}t^{1/d_w^{\mathrm{}}}`$. This claim is supported by the result $`k_ml_m`$ for $`N1`$ found by simulations of $`S_N(t)`$ and $`V(\mathrm{})`$ in two-dimensional incipient percolation aggregates, where $`l_m=\kappa _m(V)/V(\mathrm{})^m`$ characterizes the fluctuations in the volume $`V(\mathrm{})`$. Therefore, the result $`k_ml_m`$ for $`N1`$ implies that the fluctuations in $`S_N(t)`$ are mainly accounted for by the fluctuations in $`V(\mathrm{})`$, and that the fluctuations in $`S_N(t)`$ induced by the randomness of the diffusion process are irrelevant. One expects this also to be true for other disordered media. Hence, if $`S_n(t)`$ is known, the cumulants of the distribution of $`S_N(t)`$ can be calculated (for any sufficiently large value of $`N`$) from the cumulants of the distribution of the chemical volume. Finally, taking into account that $`S_N(t)`$ is reasonably well known SNtFrac , we arrive at a closed expression for the survival probability $`\mathrm{\Phi }_N(t)`$ \[c.f. Eq. (12)\] using Rosenstock’s approximation. But from $`\mathrm{\Phi }_N(t)`$ one gets the probability $`\mathrm{\Phi }_{jN}(t)`$ that $`j`$ random walkers of the initial set of $`N`$ have been absorbed by time $`t`$ \[c.f. Eq. (2)\], so that, finally, we get the moments of the trapping times, $`t_{j,N}^m`$ from the first moment $`S_n(t)`$ of the territory explored!

Comparison with simulation data shows that, in contrast with the Euclidean case, Rosenstock’s approximations of order higher than zero are necessary to account for the order statistics results in the two-dimensional percolation aggregate. This is a consequence of the large value (in comparison with the Euclidean case) of $`k_m`$, due to the large fluctuations in the territory explored by the random walkers, induced, as we showed, by the spatial disorder of the substrate. However, some features of the order statistics of trapping hold in the disordered case: for example, we found that the ratio $`\sigma _{j,N}/t_{j,N}`$ depends only on the lattice characteristic parameters $`d_s`$ and $`k_m`$, $`m=1,2,\mathrm{}`$ for large $`N`$. This was confirmed by simulations in two-dimensional percolation aggregates.

There are some interesting problems that we still cannot answer with the theory developed in this paper. For example, an important quantity is the time $`t_{N,N}`$ elapsed until all the particles of the inital set of $`N`$ are trapped. The evaluation of the moments of this quantity would require specific techniques for $`jN`$ as our results are limited to the opposite limit $`jN`$. Moreover, the trapping of the last particles surely takes place in the Donsker-Varadhan time regime DV  where Rosenstock’s approximation for the survival probability cannot be used. The recent development of a Monte Carlo method to evaluate confidently the survival probability in the Donsker-Varadhan time regime for Euclidean lattices by Barkema et al.  Barkema  following a previous work of Gallos et al.  Gallos  could serve as starting point for tackling this problem. However, one should be aware that this task is not a straightforward generalization to disordered media of that carried out for Euclidean lattices because one has to take into consideration that, as Shapir Shapir  pointed out, the Donsker-Varandhan long-time behavior is dominated by the subset of lattice realizations that are more ramified (with $`d_s=1`$). Consequently, an efficient Monte Carlo technique to explore the relevant rare lattice realizations in percolation clusters has to be devised before.

###### Acknowledgements.

This work was supported by the Ministerio de Ciencia y Tecnología (Spain) through Grant No. BFM2001-0718.


format: text

text:
# 1 Introduction

## 1 Introduction

A graph is defined as a set of nodes (vertices) and a set of links between nodes (edges) . By graph evolution or growth we mean subsequent attaching of new nodes with $`m`$ edges to previously existing nodes . Such growing graphs may reflect some features of real evolving networks, e.g. a network of collaborators, a network of citations of scientific papers, some biological networks (food chains or sexual relations) or Internet and world-wide-web pages with links between them .

The distance between nodes is the shortest number of edges which leads from one node to the other. The node-node distance (NND) distribution depends on how subsequent nodes are attached. If each node is connected with only one of preexisting nodes ($`m=1`$) a tree appears. When a new node is attached to several different nodes with $`m>1`$ edges, the growing structure is called a simple graph. We may choose nodes to which new nodes are attached in preferential way or randomly. In the latter case we deal with exponential networks. If the probability of choosing a node is proportional to its degree (e.g. to the number of its nearest neighbors) the growing structure is called scale-free or Albert–Barabási networks .

In this paper, the numerical algorithm for the network growth — basing on distance matrix evolution — is presented both, for exponential and scale-free networks ($`m=1,2`$) . The NND distribution and its characteristics are calculated. For the exponential trees the iterative formulas for $`n`$-th ordinary moments are derived.

## 2 Computer simulations

A graph with edges of unit length may be fully characterized by its distance matrix $`𝐒`$, an element $`s_{ij}`$ of which is equal to the shortest path between nodes $`i`$ and $`j`$. This matrix representation is also particularly useful when computer simulations for graph evolution are applied.

Attaching subsequent node with one edge ($`m=1`$) to previously existing network of $`N`$ nodes corresponds to adding a new $`(N+1)`$-th row and a new column to $`N\times N`$ large distance matrix $`𝐒`$. The distance from newly added $`(N+1)`$-th node to all others via selected — labeled as $`p`$ — node is larger by one than the distance between $`p`$-th and all others. Thus, new $`(N+1)`$-th row/column is a simple copy of $`p`$-th row/column but with all of its elements incremented :

$$1iN:s_{N+1}(N+1,i)=s_{N+1}(i,N+1)=s_N(p,i)+1.$$
(1)
Similarly, when new node is attached to the network with two edges ($`m=2`$) to two different nodes — which are labeled as $`p`$ and $`q`$ — but the distance from all other nodes $`i`$ to the newly added $`(N+1)`$-th is one plus the smaller distance between $`p`$$`i`$ or $`q`$$`i`$ nodes pairs :

$$1iN:s_{N+1}(N+1,i)=s_{N+1}(i,N+1)=\mathrm{min}(s_N(p,i),s_N(q,i))+1.$$
(2)
In the case mentioned above of the growth of the simple graphs also the reevaluation of distances between nodes $`i`$ and $`j`$ must be done to check if adding a new node provides the shortcut :

$$\begin{array}{c}\hfill 1i,jN:s_{N+1}(i,j)=\\ \hfill \mathrm{min}(s_N(i,j),s_N(i,p)+2+s_N(q,j),s_N(i,q)+2+s_N(p,j)).\end{array}$$
(3)
In both cases diagonal elements of new row/column are zero :

$$s_{N+1}(N+1,N+1)=0.$$
(4)
Selecting rows/columns (nodes to which we attempt to add a new node) may be random or preferential. In the latter case an additional evolving vector is introduced, which contains the node labels. These labels occur as vector elements with a probability proportional to the degree of the node. Random selection of elements of such a vector correspond to Albert–Barabási construction rule. The procedure is known as the Kertész algorithm .

## 3 Analytical calculations

Let us define $`n`$-th moments of the NND distribution for all distances

$$\mathrm{}_N^n[\{s^n(i,j)\}]=\frac{1}{N^2}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}[s^n(i,j)],$$
(5)
and only for non-zero distances

$$d_N^n[s^n(i,j)]=\frac{1}{N(N1)}\underset{i=1}{\overset{N}{}}\underset{\begin{array}{c}j=1\\ ji\end{array}}{\overset{N}{}}[s^n(i,j)],$$
(6)
where $`\{\mathrm{}\}`$, $`\mathrm{}`$ and $`[\mathrm{}]`$ denote an average over $`N^2`$ matrix elements, an average over $`N(N1)`$ non-diagonal matrix elements, and an average over $`N_{\text{run}}`$ independent realizations of the evolution process (matrices), respectively. Moments (5) and (6) for $`n=1`$ are sometimes called the network diameter. Both double sums in r.h.s. of Eqs. (5) and (6) are equal, due to obvious fact, that $`s(i,i)=0`$. That allows to derive simple dependence between averages $`\{\mathrm{}\}`$ and $`\mathrm{}`$:

$$N\mathrm{}_N^n=(N1)d_N^n.$$
(7)
For the exponential trees — basing on $`s(i,i)=0`$ and distance matrix symmetry $`s(i,j)=s(j,i)`$ — we are able to construct iterative equations for $`\mathrm{}_{N+1}^n`$ as dependent on $`\mathrm{}_N^k`$ ($`k=1,\mathrm{},n`$)

$$(N+1)^2\mathrm{}_{N+1}^n=\underset{i=1}{\overset{N+1}{}}\underset{j=1}{\overset{N+1}{}}[s^n(i,j)]=N^2\mathrm{}_N^n+2\underset{i=1}{\overset{N}{}}\left(1+[s(i,q)]\right)^n,$$
(8)
where $`q`$ is the number of the randomly selected row/column of the distance matrix $`𝐒`$. Combining Eq. (8) with Eq. (7) gives the desired iterative formula:

$$d_{N+1}^n=\frac{(N+2)(N1)}{(N+1)N}d_N^n+\frac{2}{N+1}+\frac{2(N1)}{(N+1)N}\underset{k=1}{\overset{n1}{}}\left(\genfrac{}{}{0pt}{}{n}{k}\right)d_N^k.$$
(9)
## 4 Results and conclusions

For the trees the mean of the NND $`d_N^1`$ and its dispersion $`\sigma ^2d_N^2(d_N^1)^2`$ grow logarithmically with $`N`$ (see Tabs. 1, 2) . For the graphs only the first cumulant (the average of the NND $`d_N^1`$) grows logarithmically (see Tab. 2) . Such a slow increase of $`d_N^1`$ with number of network nodes is known as the small-world effect .

The histogram of NND is presented in Fig. 1. As we expected the NND for the graphs are more condensed than the NND for the trees, as well as, the scale-free graphs (trees) are more condensed than the exponential graphs (trees).

Knowing the moments $`d_N^n`$ — the averages of the $`n`$-th powers of the non-diagonal distance matrix elements (6) — allows to build all statistical parameters which characterize the NND distribution, e.g. the average distance $`d`$, the distance dispersion $`\sigma ^2`$, its skewness

$$\nu _3\frac{d_N^33d_N^2d_N^1+2(d_N^1)^3}{\sigma ^3},$$

and kurtosis

$$\kappa _4\frac{d_N^44d_N^3d_N^1+6d_N^2(d_N^1)^23(d_N^1)^4}{\sigma ^4}.$$

The values of such characteristics of the NND for the exponential trees obtained via Eq. (9) are presented in Figs. 2 and 3. For trees, the distributions are similar to the Poisson distribution (see Fig. 1). However, even for large $`N`$ the skewness and kurtosis do not vanish as one may expect for the normal distribution .

## Acknowledgments

K.M.’s participation in the 37-th Polish Physicists’ Meeting in Gdańsk was financed by Kraków Branch of the Polish Physical Society (PTF). The numerical calculations were carried out in ACK–CYFRONET–AGH. The machine time on SGI 2800 is financed by the State Committee for Scientific Research (KBN) in Poland with grant No. KBN/SGI2800/AGH/018/2003.


text:
# BeppoSAX observations of rotation-powered pulsars

## 1 INTRODUCTION

BeppoSAX satellite operated between April 1996 and April 2001. Its scientific payload included 4 Narrow Field Instruments (hereafter NFIs; for a detailed description of the scientific payload see Boella et al. 1997) covering nearly three orders of magnitude for energy from 0.1 up to 300 keV. During its lifetime BeppoSAX observed about 20 fields of isolated pulsars detecting the source in more than half of them.

In this review, the most relevant spectral and timing results for the three more intense Crab-like pulsars, Crab, PSR B1509-58 and PSR B0540-69, and three millisecond pulsars, PSR J0218+42, PSR B1937+21 and PSR B1821-24 are reported.

## 2 THE CRAB PULSAR

The Crab pulsar (PSR B0531+21) is the only rotation powered pulsar observed at almost every energy band of the electromagnetic spectrum. Its pulse profile is characterized by a double peak structure with a phase separation of 0.4, approximately aligned in absolute phase over all wavelengths. The relative intensity, height and width of the two peaks varies with energy: in particular, the first peak (P1), dominant al low energies, becomes smaller than the second one (P2) in the soft $`\gamma `$ rays. The enhancement of the bridge between these peaks - hereafter named inter-peak (Ip) region - as function of energy is also well evident. The Crab pulsar was first detected in the X-ray band by Frits et al. (1969) and by Bradt et al. (1969). The HEAO-2 satellite produced the first high resolution ($``$ 4”) image of the Crab nebula/pulsar in X-ray along with the phase profile (Harden and Seward 1984). Pravdo et al. (1997) reported detailed pulse-phase resolved spectral analysis of the pulsed emission obtained with the PCA (5 - 60 keV) and HEXTE (16 - 250 keV) instruments on board RossiXTE. They found systematic spectral changes in the photon power law index as a function of pulse phase. Recently, Kuiper et al. (2001) by using the high energy $`\gamma `$-ray data from the CGRO satellite together with data obtained at soft/hard X-ray energies from other observatory obtained an exhaustive high-energy picture of the Crab pulsar from 0.1 keV up to 10 GeV.

BeppoSAX observed the Crab Nebula and Pulsar several times because this source was used for periodical calibration of the NFIs. All these observations provide a high statistics data set ($``$ 200 ks) that allowed to perform a very accurate timing and spectral analysis over an energy range wider than two orders of magnitude from about 0.1 up to 300 keV (see Massaro et al. 2000 for details). Timing analysis for this source was performed producing phase histograms for each NFI and each pointing. The values of P and $`\dot{P}`$ were derived from the Jodrell Bank Crab Pulsar Monthly Ephemeris (http://www.jb.man.ac.uk). The good signal to noise ratio of the NFIs, together with the high statistics allowed to describe the pulse profiles over the entire band in very high detail as shown in the top panels of Fig. 1. The morphology change of the profile as function of the energy is striking, with the second peak increasing from the soft range of the LECS (0.1-10 keV) up to the harder range of the PDS (15-300 keV).

Phase resolved spectral analysis was performed with a minimum phase resolution of 0.00667 unless the count number was not high, as in the interpeak region and in some energy ranges, where wider phase intervals were taken to reduce the statistical uncertainty in estimating the spectral parameters. Modeling the spectra with a simple power law, photon indices were computed for each NFI. A column density $`N_H`$=3.2$`\times `$1021 cm-2 derived from the spectrum of the off-pulse interval (Massaro et al. 2000) was fixed for all following fit procedures. The bottom panels in Figure 1 show the phase evolution of the spectral index in the four NFIs. The photon index changes along with the phase with a similar phase dependence in each panel: The softer emission is in P1 and P2 region while the hardest is in the middle of Ip; the spectral index difference $`\mathrm{\Delta }\alpha `$ is of the order of 0.3-0.5. Moreover, comparing the panels from left to right, we can note that the photon index, relative to the same phase interval, significantly increases with the energy. In particular, photon index of P1 changes from 1.6 in the LECS range to 2.1 at higher energies, that of P2 from 1.8 to 2.0 and Ip from 1.2 to 1.8. The phase evolution of the spectral index is very similar to that found by Pravdo et al (1997), but the spectral softening towards higher energies is not so clear as in our results because of the worse S/N ratio of the RXTE data at energies greater than 100 keV.

To explain the phase evolution of the spectral index and the change of the pulse shape with the energy a two component model was proposed by Massaro et al. (2000). They considered that the observed pulsed emission is due to the superposition of two components having different phase and energy distributions. The first component was assumed to have the same pulse profile as observed at optical frequencies with P1 much more prominent than P2 and a very low intensity in the Ip region, the second one with a harder spectrum and a phase profile estimated from the BeppoSAX data, by means of a fitting procedure of several pulse profiles at different energies, in order to reproduce the observed pulse profile when it was summed to the first component. This second component was found to be much higher than the first one in the Ip region, while the two components have a comparable intensity in the P2 region.

As a consequence of the photon index behaviour with energy, shown in Figure 1, we found that a single power law was not able to give a satisfactory representation of the spectral distribution, for a fixed phase interval region, over the entire BeppoSAX range (0.1-300 keV). Better fits were obtained using a second order law in the double log representation:

$$F(E)=K(E/E_0)^{(a+bLog(E/E_0))}$$
$`(1)`$
where $`E_0`$ is taken equal to 1 keV, and therefore $`a`$ corresponds to the photon index at this energy, while $`b`$ measures the curvature of the spectral distribution. Applying this model to three wider phase intervals, the first peak (P1) region (0.99-0.01667), the Interpeak (Ip) region (0.01667-0.28) and the second peak (P2) (0.3833-0.4167) we obtained a value for the bending parameter $`b`$ of $``$ 0.15 over the three phase intervals, and for the $`a`$ parameter the values 1.63$`\pm `$0.08 (P1), 1.31$`\pm `$0.03 (Ip) and 1.55$`\pm `$0.09 (P2). The best fit parameters are summarized in Table 1.

## 3 PSR B1509-58

PSR B1509$``$58 is one of the youngest rotation powered pulsars, with a period of about 150 ms and a period derivative of 1.5$`\times 10^{12}`$ s s-1, the highest spin-down rate of any known pulsar. It was first discovered in the soft X-rays by Einstein (Seward & Harnden 1982) and soon afterwards detected in the radio band (Manchester, Tuohy & D’Amico 1982). Later on, its pulsed emission was also detected at hard X-ray and soft $`\gamma `$-ray energies by several satellite and balloon-borne experiments: EXOSAT (0.025–11 keV; Trussoni et al. 1990), Ginga (2–60 keV; Kawai et al. 1992), SIGMA (40–300 keV; Laurent et al. 1994), Welcome-1 (94–240 keV; Gunji et al. 1994), ROSAT (0.1–2.4 keV; Greiveldinger et al. 1995), ASCA (0.5-10 keV, Saito 1998), BATSE and OSSE aboard the Compton Gamma Ray Observatory (20-5000 keV; Wilson et al. 1993, Matz et al. 1994, Ulmer et al. 1993). A detailed study of the pulse shape and spectrum in the energy range 2–200 keV with the PCA and HEXTE aboard the Rossi X-Ray Timing Explorer has been presented by Marsden et al. (1997) and Rots et al. (1998). Kuiper et al. (1999) reported definitive detection of the pulsed signal in the COMPTEL band up to 10 MeV. Moreover, they found a source compatible with the position of PSR B1509$``$58 in the skymaps of COMPTEL and EGRET between 10 and 100 MeV.

Given its high X-ray flux, second in intensity among the isolated pulsars, the BeppoSAX observation allowed to perform timing and spectral analysis at similar statistical significance as for the Crab (Cusumano et al. 2001). The pulse profile extracted by folding the MECS data (1.6-10 keV) is shown in Figure 2. The well known single pulse shape is evident. It is markedly not symmetric with the leading wing much steeper than the trailing one. No significant variation of the profile morphology has been detected inside the BeppoSAX energy band (Cusumano et al. 2001). In the spectral analysis, the pulsed and unpulsed spectra were extracted for the phase intervals 0.17-0.53 and 0.77-1.07, respectively, adopting the phase definitions of Mardsen et al. (1997). A column density of (0.91$`\pm `$0.05) $`\times 10^{22}`$ cm-2 was evaluated by analyzing the LECS and MECS spectra of the off-pulse region, assuming that the emission is mainly originating from the nebula, at the same distance as the pulsar. The X-ray and soft $`\gamma `$-ray pulsed spectra have usually been fitted with power laws, but the resulting photon indices were generally found to vary with the instrument energy windows. Trussoni et al. (1990) first measured a photon index of about 1.1 with an EXOSAT observation; the same value was found by Saito (1998) with ASCA. Ginga and RXTE observations gave a photon index of 1.35 (Kawai et al. 1993, Marsden et al. 1997), and an even steeper value of 1.68 was derived from BATSE and OSSE observations (Wilson et al. 1992, Matz et al. 1994), suggesting a softening towards the higher energies. However, no evidence for a spectral break was seen in the RXTE data up to $`200`$ keV (Mardsen et al. 1997). Figure 3 reports the photon index obtained by fitting the pulsed spectra for each of the BeppoSAX NFI with simple power laws. Notice that the photon index for the PDS is significantly larger than that for the MECS, while that derived for the LECS is smaller.

In analogy with the Crab Pulsar, a simultaneous fit of a single power law shape to the pulsed emission over the entire BeppoSAX range did not give a satisfactory result ($`\chi _\nu ^2`$ = 1.5 with 39 d.o.f.). Fitting the whole BeppoSAX spectrum of PSR B1509-52 with the continuously steepening law of Eq. (1) gives a better representation of the spectral distribution ($`\chi _\nu ^2`$ = 0.74 with 38 d.o.f.). The best fit values of the parameters were: $`a`$ 0.96$`\pm `$0.08 and $`b`$ 0.16$`\pm `$0.04 (see Table 1 for a comparison with the Crab results). The extrapolation of Eq. (1) to higher energies predicts fluxes in agreement with the values measured by OSSE, BATSE and COMPTEL (Wilson et al. 1993; Matz et al. 1993; Ulmer et al. 1993: Kuiper et al. 1999). This result is clearly shown in Figure 3. For comparison, the extrapolation of the simple power–law model is also shown.

## 4 PSR B0540-69

Another Crab-like pulsar whose pulsed emission was detected by BeppoSAX is PSR B0540-69. This pulsar, located in the Large Magellanic Cloud has a pulsed period of about 50 ms and a large period derivative of $`4.79\times 10^{13}`$ s s-1, comparable to that of the Crab pulsar. PSR B0540$``$69 was discovered in the soft X-rays by Seward et al. (1984) with the Einstein Observatory. Pulsations at optical frequencies were soon detected by Middleditch and Pennypacker (1985) with a mean pulsed magnitude of 22.5. In the radio band PSR B0540$``$69 is a quite faint source and pulsed signals were first observed only in 1989-90 (Manchester et al. 1993). BeppoSAX detected the pulsar only by the imaging LECS and MECS instruments (Mineo et al. 1999 and Cusumano et al. 2003). Several estimates of the braking index of this pulsar have been reported in the literature ranging from 2.01 $`\pm `$ 0.02 (Manchester and Peterson 1989, Nagase et al. 1990) to 2.74 $`\pm `$ 0.10 (Ögelman and Hasinger 1990). The BeppoSAX estimate of the pulsar period, combined with earlier ASCA period evaluations, allowed to measure a braking index of 2.10 $`\pm `$0.1. This results has been recently confirmed by a combined timing analysis on an extended data set including ASCA, BeppoSAX and RXTE observations spanning a time interval of about 8 years (Cusumano et al. 2003).

Given the larger distance with respect to the previous two pulsars, BeppoSAX data had not enough statistics to detected pulsed emission at energies greater than 10 keV. The energy band was not wide enough to investigate the presence of bending in the spectral distribution like that found in the other two Crab-like pulsars. Recently, de Plaa et al. (2003), making use of a long set of RXTE observations derived the pulse profiles and spectral distribution up to 50 keV. By fitting RXTE data together with the ROSAT ones they found that pulsed emission can not be described by a single power law because of the presence of a bend towards higher energies. In analogy with the results shown just above, the spectral distribution of this pulsar is well described by Eq. (1) with a value of the bending parameter similar to that found for Crab and PSR B1509-58 (Table 1).

## 5 PSR J0218+4232

PSR J0218$`+`$42 is a 2.3 millisecond pulsar in a two day orbit around a white dwarf companion (Navarro et al. 1995). The pulsar has spin a down energy of $`\dot{E}`$=2.5$`\times `$1035 erg s-1, a bipolar magnetic field component at the star surface of B=4.3$`\times `$108 G and spin down age $`\tau `$ $``$ 4.6$`\times `$108 yr. Soft X-ray emission and pulsation were first detected by ROSAT HRI (Verbunt et al. 1996; Kuiper et al. 1998) with a pulse profile characterized by a sharp main pulse and an indication of a second peak. Emission at $`\gamma `$-rays from this pulsar has been marginally detected with EGRET on board GRO with a double peak profile with $``$ 3.5 $`\sigma `$ significance (Kuiper et al. 2000). Recently, Chandra X-ray Observatory (CXO) and Rossi X-ray Timing Explorer (RXTE) observations were able to aligned the X-ray pulse profile with 2 of the 3 pulses visible at radio-frequencies and with the two $`\gamma `$-ray pulses detected in EGRET, increasing the significance of the $`\gamma `$ detection to a 4.9 $`\sigma `$ (Kuiper et al. 2002).

The BeppoSAX observation of the 2.3 millisecond pulsar PSR J0218+4232 provided for the first time detailed information on the pulsar’s temporal and spectral properties over the broad energy band 1.6–10 keV (Mineo et al. 2000). Pulsed emission was detected in the MECS energy band (1.6–10 keV) with a significance of 6.8 $`\sigma `$ at the frequency extrapolated from radio ephemeris. Figure 4 shows the light curves in the energy bins 1.6-4 keV and 4-10 keV resulting from folding all MECS events with a phase resolution of 18 bins ($``$ 0.13 ms). The pulse profiles in the MECS ranges are characterized by a double peak structure with a relative phase separation of 0.47$`\pm `$0.05. The MECS background level, determined from the spatial analysis is indicated with a dashed line.

In the same figure (top panel) the ROSAT (0.1–2.4 keV) profile (Kuiper et al. 1998) is shown, shifted in phase to obtain the highest peak coincident with the most significant one in the MECS softer light curve (middle panel).

These profiles clearly show a change of the relative peak intensities. The peak at phase 0.8 is stronger in the low energy histograms (top and middle panels), while that at phase 0.3, which is not prominent in the ROSAT low energy phase histogram, becomes the dominant feature above 5 keV (bottom panel). This behaviour reminds the P1/P2 ratio observed in the Crab light curve at higher energies. The unpulsed component determined by applying the bootstrap method (Swanepoel et al. 1996) is shown as shaded area in all three panels of Figure 4. Notice that, while in the ROSAT profile a DC component is apparent above the background level, the same does not hold for the 4–10 keV profile, where the background level is consistent with the intensity measured in the valleys of the light curve. In the intermediate energy profile there is evidence for a DC component but its intensity is not as high as in the ROSAT profile. Such effect can be interpreted as the presence of a quite soft unpulsed emission.

The pulsed emission is well modeled with a power–law, absorbed at low energy by the galactic column density with NH of 5$`\times `$1020 cm-2 (see Verbunt et al. 1996). The measured photon index is 0.61$`\pm `$0.32 harder than for any other isolated pulsar. Spectra for the pulses P1 and P2 were also fitted with power-law models with spectral indices of $`0.84\pm 0.35`$ and $`0.42\pm 0.36`$, respectively, in agreement with the trend seen in Figure 5. The unabsorbed (2–10 keV) pulsed flux is 4.1$`\times `$10-13 erg cm-2 s-1 implying a luminosity of $`L_x=1.3\times 10^{32}\mathrm{\Theta }(d/5.7\mathrm{kpc})^2`$ erg s-1, where $`\mathrm{\Theta }`$ is the solid angle spanned by the emission beam.

## 6 PSR B1937+21

PSR B1937$`+`$21, the fastest known millisecond pulsar (P=1.558 ms), is located at a distance of 3.6 kpc as based on dispersion measure (Taylor & Cordes 1993). Its spin down energy is $`\dot{E}`$=1.1$`\times `$1036 erg s-1 and the dipolar magnetic field component at the star surface is B=4.1$`\times `$108 G.
It was discovered in radio by Backer et al. (1982). X–ray emission from this pulsar was first detected by ASCA (Takahashi et al. 2001) above 2 keV, with a pulse profile characterized by a single sharply peak and a pulsed fraction of 44%. ASCA data allowed to aligned the X-ray peak to the radio interpeak.

BeppoSAX detected pulsed emission from PSR B1937+21 (Nicastro et al. 2002), the fastest known millisecond pulsar (1.55 ms). Pulsed emission was detected in the MECS energy (1.6–10 keV) with a significance of 11 $`\sigma `$ at a frequency deviating from that extrapolated from the radio ephemeris by $`7.6\times 10^6`$ Hz. This discrepancy is consistent with a systematic frequency error characteric for the timing precision of the BeppoSAX clock for observation later than 2000 (Mineo et al. 2003). Figure 5 shows the X-ray pulse profiles for the whole MECS range (top panel) and for two energy sub-intervals 1.6–4 (middle panel) and 4–10 keV (bottom panel). The pulse profile is characterized by a double peak with a phase separation of $`0.48\pm 0.04`$. The significance of the second peak, first detected by BeppoSAX, is about 5 $`\sigma `$. The comparison of the middle and bottom panels in Figure 5, indicates a hint for variations in the ratio of the two peaks. The unpulsed level determined by applying the bootstrap method and the background level derived from a spatial analysis are also shown in Figure 6. The DC level is compatible with the background level implying that the source is consistent with being 100% pulsed.

The pulsed emission was modeled by an absorbed power law with a photon index of 1.66$`\pm `$0.07 and $`N_\mathrm{H}`$ of $`2.3\times 10^{22}`$ cm-2. The absorbed flux (2–10 keV) is $`3.1\times 10^{13}`$ erg cm-2 s-1 implying a luminosity of $`L_X=4.6\times 10^{31}\mathrm{\Theta }`$ ($`d/3.6`$ kpc)2 erg s-1 and an X-ray efficiency of $`\eta =4\times 10^5\mathrm{\Theta }`$.

## 7 PSR B1821-24

PSR B1821$``$24 is an isolated millisecond radio pulsar discovered in the core of the globular cluster M28 (Lyne et al. 1987) with high spin-down luminosity ($`\dot{E}`$=2.2$`\times `$1036 erg s-1). It has a spin down period of P=3.05 ms and the dipolar magnetic field component at the star surface is B=2.2$`\times `$109 G. A marginal detection of X-ray pulsation was first reported by Danner, Kulkarni and Thorsett (1994) from a ROSAT (PSPC) observation. A more significant detection was performed by ASCA (Saito et al. 1997). The pulse profile is characterized by a double peak with a phase separation of about 0.44. RXTE and Chandra observations allowed to aligned the narrower X-ray peak to one of the radio peak (Rots et al. 1998; Rutledge et al. 2003). Moreover, Chandra data provided detailed spectral results on the pulsed emission (Becker et al. 2002).

BeppoSAX observation of PSR B1821-24 has provided more accurate information on the timing and spectral properties of the pulsar. The pulsation in the MECS range (1.6 – 10 keV) was detected at high significance level (15 $`\sigma `$) at a frequency in agreement with the radio ephemeris within the systematic deviation of the onboard clock (Mineo et al. 2003). Figure 6 shows the X-ray pulse profile in the total MECS range (top panel) and in two energy sub-intervals 1.6–4 (middle panel) and 4–10 keV (bottom panel). The shaded areas indicate the DC level ($`\pm 1\sigma `$), while the dashed lines indicate the background level. Note that the DC level shown in the figure includes the contribution of the M28 globular cluster in which the pulsar is located. Also for this millisecond the pulse profile is characterized by a double peak structure with phase separation of 0.45$`\pm `$0.03. A hint for variation in the ratio of the two peaks with energy can be inferred comparing the middle panel with the bottom panel of Figure 6.

The pulsed emission is well modeled with a power–law, absorbed at low energy by the galactic column density with NH of 2.8$`\times `$1020 cm-2 and spectral index of 1.45$`\pm `$0.21. The unabsorbed (2–10 keV) pulsed X-ray flux is 3.9$`\times `$10-13 erg cm-2 s-1 implying a luminosity of $`L_x=1.4\times 10^{33}\mathrm{\Theta }(d/5.5\mathrm{kpc})^2`$ erg s-1.

## 8 CONCLUSION

On the basis of high quality BeppoSAX data, covering nearly three orders of magnitude in energy (0.1-300 keV) we have showed that the energy spectra of the most intense Crab-like pulsars can not be described by a single power lawn shape, but they show a significant steepening towards higher energies. In fact, modeling the spectra of Crab, finely selected in phase, with a simple power law we obtained photon indices significantly increasing with the energy window of the NFIs. A similar result was also obtained for PSR B1509-58. Broad band spectra of three phase intervals of the Crab light curve (P1,P2 and IP) and of the pulsed emission of PSR B1509-58 can be well modeled by a curved power law, Eq. (1), characterized by a continuously bending of the spectral energy distribution. Recently, De Plaa et al. (2003) found, analyzing a set of RXTE observations, that PSR B0540-69 shows a similar spectral behaviour. The $`b`$ parameter, representing the curvature of the spectral distribution, has been found to be similar, equal $``$ 0.15 for all spectra. We stress that this new result has been obtained just thanks to the broad band capability of BeppoSAX and the excellent inter-calibration among the 4 Instruments.

Such curved spectra suggest the existence of a possible characteristic energy, for instance the peak of the Spectral Energy Distribution (SED). In the scenario of the outer gap model (Cheng et al. 1986a, 1986b) the energy of the maximum of the SED could be easily related, at least for this three Crab-like pulsars, by a simple scaling law to some observable pulsar parameters as the period and its first derivative (Cusumano et al. 2001):

$$(E_m/E_{mC})(P/P_C)^{69/40}(\dot{P}/\dot{P}_C)^{11/40}$$

$$(P/P_C)^{7/4}(\dot{P}/\dot{P}_C)^{1/4}.$$
$`(2)`$
where the index C refers to the Crab values. The maximum of the SED has been calculated for the Crab to be 14 keV and 200 keV for the P1 and Ip region, respectively, for PSRB1509-58 a higher value around 5 MeV (Cusumano et al 2001), and for PSR B0540-69 the maximum is located around 170 keV.

The above scaling relation of the energy of the maximum of the SED of PSR B1509-58 and PSR B0540-69 with respect the maximum for the Crab matches quite well when the spectrum of the Crab interpeak region is considered, while it is about one order of magnitude lower for that of the first peak of Crab. This suggests a similar mechanism of energy production for the X-ray emission of the PSR B1509-58, PSR B0540-69 and for the Crab interpeak component that correspond to the second component of the light curve in the model proposed by Massaro et al. (2000).

Pulsed emission has been detected in the range 1.6-10 keV for the 3 millisecond pulsars. All detected light curves have double peak profiles with high pulsed fraction and with a phase separation of about 0.5. PSR J0218+4232 shows spectral variation with phase: the relative intensity of the two peaks varies with energy with a behaviour that strongly recalls the Crab pulsar. Similar behaviour seems also present in the other two milliseconds (PSR B1937+21 and PSR B1821-24), although it needs to be confirmed with future more significant data set. The pulsed emission of these milliseconds is modeled by a flat power law indicating a likely non-thermal emission. The X-ray profile morphology and the hard non-thermal spectra suggest a similar magnetospheric origin of the high-energy emission of these millisecond pulsars and the Crab. The comparison with Crab is additionally enforced by the similar values of the magnetic field at the light cylinder as already pointed out by Kuiper et al. (1998) and Kawai and Saito (1999): the magnetic field strength $`B_L`$ of PSR J0218+4232, PSR B1821-243 and PSR B1937+21 is in the range of 3 and 10 $`\times `$ 105 Gauss, and that of the Crab is 9 $`\times `$ 105 Gauss. On the other hand, the magnetic field strength at the neutron star surface $`B_S`$ of the millisecond pulsars is more than 4 order of magnitude weaker than the Crab one.

The similar value of $`B_L`$ seen in these millisecond pulsars and in the Crab suggest that the magnetic field strength near the light cylinder is a key parameter to explain their high energy emission, indicating as likely model for the high energy emission the outer gap model. However, the polar cap scenario can also account for high-energy emission from these milliseconds pulsars (Luo, Shibata and Melrose 2000; Dyks and Rudak 2002).

## ACKNOWLEDGEMENTS

I am grateful to T. Mineo, E. Massaro and L. Nicastro for their useful comments and their help for preparing this paper. The author is very grateful to the referee, W. Hermsen, for the comments and suggestions that improved the paper.

E-mail adress of G. Cusumano cusumano@pa.iasf.cnr.it


text:
# Comment on “Entanglement Measure for Composite Systems”

## Abstract

A recent paper claimed to give an entanglement measure for composite systems, including Bose condensation and superconductivity. It is not an entanglement measure. It does not distinguish entanglement and non-factorization merely due to (anti)symmetrization for identical particles. In fact, the nature of entanglement, beyond the effect of (anti)symmetrization, in important states of many identical particles such as Bose condensation and superconductivity, had already been studied earlier.

A recent paper yukalov  claimed to give an entanglement measure for composite systems. It came to our notice that what was defined is not an entanglement measure. Moreover, for identical particles, this measure does not distinguish entanglement and non-factorization merely due to (anti)symmetrization, e.g. the Hatree-Fock states, which is not entangled shi2 .

In yukalov , a measure $`ϵ(A)`$ was defined for an arbitrary operator $`A`$ acting on a composite system, which vanishes when $`A`$ is a tensor product of operators on subsystems. Then the operator was taken to be the density matrix $`\rho `$, and $`ϵ(\rho )`$ was used to measure whether $`\rho `$ is entangled. Several examples, including Bose condensation and superconductivity, were discussed by using this measure.

It can be seen that while $`ϵ(\rho )`$ is indeed nonzero for an entangled pure state of distinguishable particles, it is nonzero also for an entirely classical mixture.

Moreover, $`ϵ(\rho )`$ is nonzero for symmetrized or antisymmetrized non-entangled state of a system of identical particles when they occupy more than one single particle states. In fact, entanglement in a system of identical particles should be beyond the effect of (anti)symmetrization peres ; shi1 ; zanardi .

In none of these two cases, the nonzero $`ϵ(\rho )`$ measures the entanglement. This is contrary to the idea in yukalov .

Furthermore, the nature of entanglement, beyond the effect of (anti)symmetrization, in important states of many identical particles, such as Bose condensation and superconductivity, had already been studied in detail earlier shi2 ; shi3 .


text:
# Partially incoherent optical vortices in self-focusing nonlinear media

## Abstract

We observe stable propagation of spatially localized single- and double-charge optical vortices in a self-focusing nonlinear medium. The vortices are created by self-trapping of partially incoherent light carrying a phase dislocation, and they are stabilized when the spatial incoherence of light exceeds a certain threshold. We confirm the vortex stabilization effect by numerical simulations and also show that the similar mechanism of stabilization applies to higher-order vortices.

Vortices are fundamental localized objects which appear in many branches of physics, ranging from liquid crystals to superfluids and Bose-Einstein condensates pismen . In optics, vortices are associated with phase dislocations (or phase singularities) carried by diffracting optical beams soskin . When such vortices propagate in self-defocusing nonlinear media, the vortex core with a phase dislocation becomes self-trapped, and the resulting stationary singular beam is known as an optical vortex soliton book . Such vortex solitons have been generated experimentally (e.g. by using a phase mask) within broad diffracting beams in different types of defocusing nonlinear media swartz ; exp1 ; exp2 ; exp3 ; review . Optical vortex solitons demonstrate many common properties with the vortices observed in superfluids and, more recently, in Bose-Einstein condensates bec .

In contrast, optical vortices become highly unstable in self-focusing nonlinear media. Indeed, when a ring-like optical beam with zero intensity at the center and a phase singularity kruglov  propagates in a self-focusing nonlinear medium, it always decays into several fundamental solitons flying off the main ring firth . This effect has been observed experimentally in different nonlinear systems, including the saturable Kerr-like nonlinear media tikh , biased photorefractive crystals exp3 , and quadratic nonlinear media chi2  in the self-focusing regime. This effect has also been predicted to occur in many other physical systems including attractive Bose-Einstein condensates saito .

A number of recent theoretical papers cubic , including the rigorous studies of the vortex stability buryak , suggest that the stable propagation of spatial and spatiotemporal vortex-like stationary structures may become possible in the models with competing nonlinearities in the presence of a large higher-order defocusing nonlinearity, but no realistic physical systems to support these theoretical findings have been found so far. Thus, the main question remains open: Can stable optical vortices readily be observed in self-focusing nonlinear media?

The main purpose of this Letter is to demonstrate, for the first time to our knowledge, that stable propagation of spatially localized optical vortices in a self-focusing photorefractive crystal can indeed be observed provided such vortices are created by partially incoherent light carrying a phase dislocation. In particular, we show, both experimentally and theoretically, that single- and double-charge optical vortices can be stabilized in self-focusing nonlinear media when the value of the spatial incoherence of light exceeds a certain threshold, and these vortices are readily observed in experiment as stationary self-trapped structures propagating for many diffraction lengths.

We should mention that the generation and properties of singular optical beams created by partially incoherent light is an important issue which is a subject of a current active research even in linear optics (see, e.g, Ref. gbur  and references therein). Our results demonstrate that partially incoherent singular beams can exist in self-focusing nonlinear media being stabilized by the light incoherence effect.

First, we present our experimental results. The main purpose of our experiments conducted for a biased photorefractive medium is to generate partially incoherent vortices and vortex solitons, and then inspect their stability in such a self-focusing nonlinear medium. Being driven by the earlier results of Anastassiou et al. inst\_exp  who demonstrated the suppression of modulational instability for the stripe spatial solitons created by partially incoherent light, we try to understand whether the other type of the nonlinearity-driven instability, the so-called azimuthal instability, can be suppressed by reducing the coherent properties of optical vortex beams.

The experiment (the setup is shown in Fig. 1) is conducted with a biased photorefractive SBN crystal (a x b x c = 5 mm x 10mm x 5mm). First, a cw laser light beam (at 488 nm) of the extraordinary polarization is made partially incoherent by passing it through a lens and then through a rotating diffuser. The rotating diffuser introduces random-varying phase and amplitude on the light beam every 1 $`\mu `$s, which is much shorter than the response time (about 1 s) of the crystal under our experimental light illumination. By adjusting the position of the diffuser to near (away from) the focal point of the lens in front the diffuser, we can increase (decrease) the degree of the light coherence. We collect the light after the rotating diffuser by a second lens and then pass is through a computer-generated hologram to imprint a vortex phase (with a single-, for $`m=1`$, or double-, for $`m=2`$, charge) on the light beam. Since the partially incoherent light beam can be considered as a superposition of many mutually-incoherent light beams, the first-order diffracted light beam after the hologram becomes a superposition of many mutually-incoherent vortex beams. We then focus and launch the partially coherent vortex beam into the SBN crystal along its a-axis. The crystal is illuminated by a halogen lamp from its top side as necessary for photorefractive screening solitons. The total power of the vortex beam is of 0.17$`\mu `$W, which results in the nonlinearity of the photorefractive crystal falling into the Kerr region for the peak intensity of the vortex beam to the background intensity is much less than unity. Then, a lens is used to project the images at the input and output faces onto a CCD camera.

Before showing that reducing the degree of coherence can stabilize the vortex beam propagating in a self-focusing medium, we reproduce the experiment that the coherent single-charge($`m=1`$) vortex light beam cannot stably propagate in a self-focusing medium tikh ; exp3 . We first remove the diffuser from the experimental setup. The vortex beam at the input face of the crystal is shown as Fig. 2(a). With zero biasing voltage, Fig. 2(b) shows the natural diffraction of the vortex light beam. While a 2.5 kV biasing voltage is applied on the photorefractive crystal creating a Kerr-type self-focusing nonlinear medium, the vortex beam breaks up into two pieces \[Fig. 2(c)\]. This vortex break-up observed in a self-focusing medium is due to the azimuthal instability, and it has been theoretically and experimentally demonstrated previously firth ; tikh . We then put back the rotating diffuser and adjust it to a suitable position. The degree of coherence of the vortex light beam can be estimated by the speckle size at the input face of the crystal \[Fig. 2(d)\] when we stop the diffuser from rotating. As a voltage of 2.5 kV is applied on the crystal, Fig. 2(g) clearly shows that the vortex light beam is stabilized by the reduction of the degree of coherence though two very unclear bright spots still can be seen on the opposite sides (top and bottom) of the ring-like intensity distribution. As we move the rotating diffuser further away from the focal point of the lens to make the light more incoherent \[indicated by Fig. 2(h)\] and apply a voltage of 2.5 kV on the crystal, Fig. 2(k) shows the generated stable partially incoherent vortex soliton at the output face of the crystal.

We notice that the vortex stabilization by the light incoherence can possibly be explained by employing simple physics. Indeed, after some propagation of the partially incoherent beam with an imprinted phase, we observe that the intensity at the center of the beam does not vanish. This means that the incoherent vortex can be decomposed into many mutually incoherent vortex beams with not only their phases randomly distributed but also their central positions offset from each others. In this way, the core of a composite vortex beam will be filled by some light, and the index change at the beam center will become nonzero. We believe that this effect contributes strongly to the vortex stabilization.

We continue with the experiments for the double-charge vortex beams generated by a phase mask with the $`m=2`$ dislocation. Figures 3(a-c) show that a coherent double-charge vortex light beam cannot stably propagate in a self-focusing medium (we applied 3 kV here) as well, and it breaks up into pieces as been observed in a self-focusing atomic vapor tikh . When we make the light more incoherent, the double-charge vortex becomes more stable, as shown in Figs. 3(g) and (k). However, we could not make a double-charge vortex soliton here. The double-charge vortex beam diffracts more than a single-charge vortex beam and, therefore, it requires higher nonlinearity (or higher voltage) to form the vortex soliton, but our crystal begins to arc at 3.5 kV. Nevertheless, these two experiments show that the reduction of the degree of coherence of the light beam can indeed stabilize the single- or double-charge (or even higher-charge) vortex solitons propagating in a self-focusing medium.

We have also studied numerically the propagation of partially incoherent optical vortices in a photorefractive medium, employing the coherent density approach IncNum ; inst\_exp . The coherent density approach is based on the fact that partially incoherent light can be described by a superposition of mutually incoherent light beams that are tilted with respect to the $`z`$-axis at different angles. One thus makes the ansatz that the partially incoherent light consists of many coherent, but mutually incoherent light beams $`E_j`$: $`I=_j|E_j|^2`$. By setting $`|E_j|^2=G(j\vartheta )I`$, where

$$G(\theta )=\frac{1}{\sqrt{\pi }\theta _0}e^{\theta ^2/\theta _0^2}$$
(1)
is the angular power spectrum, one obtains a partially incoherent light beam whose coherence is determined by the parameter $`\theta _0`$, i.e. less coherence means bigger $`\theta _0`$. Here, $`j\vartheta `$ is the angle at which the $`j`$-th beam is tilted with respect to the $`z`$-axis. For our simulations we used a set of 1681 mutually incoherent vortices, all initially tilted at different angles. To simulate the photorefractive nonlinearity we use a simple model which predicts that the refractive index change is approximately given by $`I/(1+I)`$.

Figure 4 shows a comparison between the propagation of the vortices generated by coherent and partially incoherent light (for $`\theta _0=0.38`$). Increasing the incoherence further leads to the case where the vortex beam radiates off a lot of its intensity before the azimuthal instability can set in. Larger values of incoherence (i.e. larger values of $`\theta _0`$) correspond to a complete suppression of the vortex instability, and this confirms the experimental results presented above.

Figures 5(a-d) show our numerical results for the propagation of an input Gaussian beam carrying a phase dislocation \[(a)\] after the total propagation (9 mm) in a nonlinear medium for the coherent light \[(b)\] and two different partially incoherent beams \[(c,d)\], corresponding to the values $`\theta _0=0.14`$ and $`\theta _0=0.29`$, respectively. The most obvious difference to the scenario of the propagation the a coherent vortex is that the vortex decay undergoes a visible delay when the degree of incoherence grows. Furthermore, in the incoherent case the vortex changes its profile only very slowly as it propagates and thus can be considered as being in a transition stage between the decay and stabilization.

Finally, we have studied numerically the propagation of a double-charge vortex beam in both coherent and partially incoherent cases, as shown in Fig. 6. The initial perturbation in both cases have been chosen to be very small in order to obtain a clearer picture of the instability, therefore the propagation distances are quite long. In numerics, since we have a possibility to observe long propagation distances, we did not find that the vortices can be completely stabilized by increasing incoherence. When using the values of $`\theta _0`$ far above the value $`\theta _0=0.35`$ used in Fig. 6, we observe that the vortex just radiates off a lot of its intensity and then decays. The problem is obviously that the vortices are not only incoherent along the angular, but also the radial direction. We believe that the experiment would confirm this observation if a longer crystal was available.

In conclusion, for the first time to our knowledge, we have observed the stable propagation of single- and double-charge optical vortices in a self-focusing nonlinear medium. The vortices have been created by partially incoherent light beams carrying a phase dislocation and propagating in a photorefractive nonlinear medium. The vortex azimuthal instability in a self-focusing nonlinear medium was suppressed for the light incoherence above a critical value. The experimental results have been confirmed by numerical simulations which also provide an insight into the physical mechanisms of the vortex stabilization.

Yuri Kivshar thanks the Physics Department of the Taiwan University for hospitality. This work was supported by a collaboration program between the Australian Academy of Science and the National Science Council of Taiwan, and the Australian Research Council.