text:
# Blue horizontal branch stars in the Sloan Digital Sky Survey: II. Kinematics of the Galactic halo

## 1 Introduction

Among all visible components of the Galaxy, the stellar halo is the oldest and occupies the largest volume. It is therefore of interest both as a record of the early history of the Galaxy and as a tracer of the Galactic potential.

The kinematics of halo stars in the solar neighborhood are well established, as shown in Table 1. For present purposes, the solar neighborhood extends as far as a few kpc from the Sun, but to a distance less than the Sun’s own Galactocentric radius. Evidently, the systemic rotation $`v_{\mathrm{halo}}`$ of the local halo is small or zero, and the velocity ellipsoid is distinctly elongated in the radial direction.

The more distant halo has been less well studied. Proper motions are generally not measurable for distant stars, so that line-of-sight velocities must be analyzed under *a priori* assumptions about the symmetries of the velocity ellipsoid. The stars are also fainter, reducing the accuracy of the measured quantities.

Much work on the distant halo has been done using globular clusters, which have the advantage of being distinctive and bright, but have important drawbacks: there are only $`10^2`$ of them, so that statistical uncertainties are large, and special conditions of formation or selective destruction by Galactic tides may have biased the properties of the cluster system compared to those of the field. Frenk & White (1980) analyzed the kinematics of a sample of 66 globular clusters under the assumption that the systemic rotation $`v_{\mathrm{halo}}`$ is constant; they concluded that the residual velocities are consistent with being isotropic within broad uncertainties. For a subsample of 44 metal-poor (“F”) clusters, they found $`v_{\mathrm{halo}}=42\pm 33\mathrm{km}\mathrm{s}^1`$ and a root-mean-square (rms) residual line-of-sight velocity $`\sigma _{\mathrm{los}}=124\mathrm{km}\mathrm{s}^1`$. In a similar analysis of clusters with $`\text{[Fe/H]}<0.8`$, Armandroff (1989) found $`v_{\mathrm{halo}}=46\pm 29\mathrm{km}\mathrm{s}^1`$ and $`\sigma _{\mathrm{los}}=116\pm 10\mathrm{km}\mathrm{s}^1`$. The restriction to low metallicities follows the discovery by Zinn (1985) that the more metal-rich clusters have a distribution and kinematics intermediate between those of the stellar halo and thin disk.

Ratnatunga & Freeman (1989) sampled the halo field using relatively faint red giants at high Galactic latitude. Their summary for the more metal-poor stars near the south Galactic pole concludes that $`\sigma _z`$ is approximately constant at $`75\mathrm{km}\mathrm{s}^1`$ out to $`25\mathrm{kpc}`$, which would be consistent with the solar neighborhood (Table 1) if the velocity ellipsoid is constant in cylindrical coordinates; however, for their most distant subsample near the SGP ($`D=17.1\mathrm{kpc}`$), they measured $`\sigma _{\mathrm{los}}=42\pm 10`$ (10 stars).

In a study similar to that reported here, Sommer-Larsen et al. (1997) analyzed a sample of 679 blue horizontal branch stars drawn from multiple sources and spanning Galactocentric radii $`765\mathrm{kpc}`$. With line-of-sight velocities and photometric distances, they found that the radial velocity dispersion (in spherical coordinates) falls rather sharply from its value in the solar neighborhood, $`\sigma _r=140\pm 10\mathrm{km}\mathrm{s}^1`$, to an asymptotic value $`89\pm 19\mathrm{km}\mathrm{s}^1`$ at $`r20\mathrm{kpc}`$. By requiring dynamical equilibrium in an assumed logarithmic Galactic potential, they concluded that the tangential velocity dispersion should rise from its local value to $`\sigma _\theta \sigma _\varphi 137\mathrm{km}\mathrm{s}^1`$ at large $`r`$. They did not claim, however, to have detected this rise directly in their data. They further suggested that the shift from radial to tangential anisotropy with increasing Galactocentric radius might support theories in which the halo formed by heirarchical merging of smaller systems.

Because of the differences in kinematics cited above, it appears that local high-velocity, low-metallicity stars may not be entirely representative of the more distant halo. To study the more distant halo, we have isolated a sample of $`1170`$ blue horizontal branch (BHB) stars from the Sloan Digital Sky Survey (Sirko et al., 2003, hereafter Paper I). In that paper, we estimated contamination from foreground stars and blue stragglers to be less than $`10\%`$ for bright stars ($`g<18`$) and about $`25\%`$ for faint stars ($`g>18`$). An improved color-absolute magnitude relation was also derived, and these absolute magnitudes are used throughout this paper to derive distance moduli. This sample of BHB stars provides effective tracers of the Galactic halo, with no kinematic bias. The sample may have selection bias, and is probably far from complete, but these biases do not affect the kinematic analyses presented in this paper other than degrading the statistics and increasing the errors. Using maximum-likelihood methods, we fit for the mean rotation and principle axes of the halo velocity ellipsoid on the assumption that these are constant in spherical or cylindrical coordinates. We also fit for the orbital velocity of the Sun.

The present paper is organized as follows. We first discuss an overview of the kinematic analysis in §2 and appendix A. The results for the kinematic parameters of the halo are presented in §3 and summarized in Table 2. Monte Carlo simulations which validate the analysis are presented in §4. §5 summarizes our main conclusions and discusses the agenda for future work.

### 1.1 Terminology used in this paper

The following terminology and conventions are used throughout this work. A heliocentric right-handed Cartesian coordinate system $`(x,y,z)`$ is oriented so that $`+x`$ points towards the Galactic center (GC) and $`+z`$ points towards the North Galactic Pole, such that $`+y`$ is the direction of disk rotation. Thus, $`𝐯_{}=(v_{,x},v_{,y},v_{,z})=𝐯_{\mathrm{LSR}}+(U,V,W)`$ where $`(U,V,W)`$ is the solar motion with respect to the local standard of rest (LSR). The velocity of the LSR is taken to be $`𝐯_{\mathrm{LSR}}=(0,220,0)\mathrm{km}\mathrm{s}^1`$ (Kerr & Lynden-Bell, 1986), and the canonical solar motion vector is $`(U,V,W)=(10,5,7)\mathrm{km}\mathrm{s}^1`$ (Dehnen & Binney, 1998, rounded to integer values). Galactocentric spherical $`(r,\theta ,\varphi )`$ and cylindrical $`(R,\varphi ,z)`$ coordinate systems will also be used. We make a rudimentary model of halo bulk rotation parametrized by a single velocity $`v_{\mathrm{halo}}`$. A positive $`v_{\mathrm{halo}}`$ indicates a flat halo rotation curve in the same sense as the disk rotation. Radial velocities to stars are measured as redshifts which are subsequently and immediately corrected for the Earth’s motion around the Sun and the Earth’s rotation about its own axis. This heliocentric radial velocity $`v_r`$ contains a component due to the solar velocity $`𝐯_{}`$.111Here the subscript “$`r`$” is used for components along the line of sight, but elsewhere it will be used for the radial component in spherical coordinates based on the Galactic Center. The context should make our meaning clear. The “line-of-sight” velocity $`v_{\mathrm{los}}`$ will be defined as $`v_r`$ plus the projection of $`𝐯_{}`$ onto the line-of-sight to the star, minus the projection of its expected rotation component $`v_{\mathrm{halo}}\widehat{\mathit{\varphi }}`$ (if any). That is, $`v_{\mathrm{los}}`$ is the component of the purely random component of velocity in a model of the halo with a flat rotation curve, projected onto the line of sight. We find $`v_{\mathrm{los}}`$ conceptually useful, but one should not forget that it is contingent upon our assumptions concerning the rotation of the Sun and the halo, whereas $`v_r`$ is directly observed (modulo errors). Throughout, the distance from the Sun to the GC is assumed to be $`8\mathrm{kpc}`$ (Reid, 1993).

## 2 The velocity ellipsoid of the Galactic halo: overview

The velocity ellipsoid is defined such that its axes align with the principal axes of the velocity dispersion tensor, and the lengths of its axes $`𝝈`$ are equal to the components of the velocity dispersion along those directions. It is difficult to know the orientation of the velocity ellipsoid so it is reasonable to assume it is oriented along galactocentric spherical or cylindrical coordinates at every point in space (Binney & Merrifield, 1998). Even with $`10^3`$ stars, it is difficult to assess whether and how the components of the velocity ellipsoid vary in space, so they are assumed constant, although we test this below by dividing our sample into a few magnitude bins.

A maximum likelihood technique is used, as outlined in Appendix A. A uniform prior is adopted (Lupton, 1993): when the probability of observing the set of stars at the observed positions and heliocentric radial velocities given an assumed $`𝝈`$ is maximized as a function of $`𝝈`$, the probability of $`𝝈`$ being the true velocity ellipsoid is maximized. In symbols, $`\varphi (v_{\mathrm{los}}|𝝈)\varphi (𝝈|v_{\mathrm{los}})`$ where $`\varphi `$ represents probability density. Naturally, $`\varphi (v_{\mathrm{los}}|𝝈)`$ depends implicitly upon the orientation of the velocity ellipsoid with respect to the line of sight, and hence upon the angular position $`(\alpha ,\delta )`$ and distance of each star. It is straightfoward to expand the maximum likelihood procedure to include four additional parameters: the three components of the velocity of the Sun with respect to the galactic center, $`𝐯_{}`$, and the bulk rotation velocity of the halo, $`v_{\mathrm{halo}}`$.

We assume that the bulk rotation of the halo corresponds to a uniform velocity in the direction of galactic rotation. It is difficult to apply this analysis on this sample to more complex halo rotation laws because of the moderate degeneracy between the halo rotation and $`v_{,y}`$ (see figure 2, discussed below); this degeneracy would increase if one attempted to fit for more parameters of the bulk halo rotation.

Distance moduli are derived from the color-absolute magnitude relation of Paper I. The measured angular coordinates of the stars $`(\alpha ,\delta )`$ are virtually error-free. However, the magnitudes, and thus distances, are susceptible to several sources of error: extinction, photon statistics, and intrinsic brightness variations from star to star which create deviations in the stars’ distance moduli. Therefore, magnitude errors should be included in the likelihood analysis. As shown in Paper I, typical magnitude errors are about $`\mathrm{\Delta }m0.2`$. In anticipation of our finding, discussed below in §3, that the kinematic analysis is very robust with respect to magnitude errors, we adopt an even more conservative magnitude error: $`\mathrm{\Delta }m=0.5`$, and assume the errors follow a gaussian distribution. Equation A3 shows the probability density $`\varphi `$ of observing a star at a given position in space, line of sight velocity $`v_{\mathrm{los}}`$, and $`𝝈`$. Having assumed that magnitude errors have a gaussian distribution with some width $`\mathrm{\Delta }m`$, the likelihood of observing a star considering magnitude errors can be found by integrating the probability density of equation A3, multiplied by a gaussian weight function for the magnitude error, along the line of sight. Figure 1 shows an illustration of this procedure for two stars with artificial data. These two stars are constructed to have $`(l,b)=(0,30)`$, $`g=15.5`$, and absolute magnitude $`0.7`$, which places them approximately 9 kpc from the Sun above the galactic center. The evaluation of the probability density $`\varphi `$ for each star alone along the line of sight is shown in the upper right panel. One star’s line-of-sight velocity is $`v_{\mathrm{los}}=0\mathrm{km}\mathrm{s}^1`$ (solid line and filled dots) and the other’s is $`v_{\mathrm{los}}=100\mathrm{km}\mathrm{s}^1`$ (dotted line and open circles). We apply $`𝝈=(150,50,100)\mathrm{km}\mathrm{s}^1`$ (spherical coordinates) to this model (note that the $`\varphi `$ component of the velocity ellipsoid is irrelevant for this particular case). It is easy to see several trends in the upper right plot of figure 1. If there were no distance measurement, for example, the stationary star would be most likely to be found at $`7\mathrm{kpc}`$ because the line of sight coincides there with the smallest ($`\theta `$) axis of the velocity ellipsoid. In fact, the probability density of observing a stationary star at that point is just $`(2\pi \sigma _\theta ^2)^{1/2}.008\mathrm{s}\mathrm{km}^1`$, as shown. Also, the receding star is least likely to be observed there because its $`v_{\mathrm{los}}`$ can most easily be “explained” by most of its contribution coming from the radial component of the velocity ellipsoid. Also note that the stationary star always has a higher probability than the receding star; this follows from the gaussian nature of the probability of observing $`v_{\mathrm{los}}`$ in equation A3.

To incorporate the magnitude error $`\mathrm{\Delta }m`$ into the analysis, the probability density $`\varphi `$ as a function of distance must be multiplied by a gaussian centered on the measured magnitude with width $`\mathrm{\Delta }m`$, and then integrated over distance. This integration is most easily accomplished with Gauss-Hermite integration, which samples $`\varphi `$ at $`N_\mathrm{p}=11`$ points, shown in the upper right panel of figure 1 as dots and circles. These 11 values of $`\varphi `$ at the specified abscissae are multiplied by the coresponding weight factors $`w_j`$ (Abramowitz & Stegun, 1972), which are plotted in a discrete fashion as the lower line in the middle panel of figure 1. The consideration of magnitude errors also introduces Malmquist bias. If the number density of observable BHB stars in the Galactic halo follows a spherical power law $`\rho r^\alpha `$ in Galactocentric radius, then Malmquist factor $`M(r^{})r^{\mathrm{\hspace{0.17em}3}}r^\alpha `$ must weight the probability density $`\varphi `$, in addition to the gaussian weight $`w_j`$ due to the magnitude error, where $`r^{}`$ is heliocentric distance. These weighting factors are normalized so that $`_j^{N_\mathrm{p}}w_jM(r)=1`$. The function $`M(r^{})`$ is shown in the middle panel of figure 1 as the top line, assuming $`\alpha =3.5`$. Note that for $`\alpha =3`$, $`M(r8\mathrm{kpc})\mathrm{constant}`$.

Finally, the probability density of observing a star at a given angular position, distance modulus, and line-of-sight velocity, with assumptions on $`𝝈`$, $`\mathrm{\Delta }m`$ and $`\alpha `$ for the Malmquist term, is just $`_j\varphi (r)w_jM(r)`$. The term $`\varphi (r)w_jM(r)`$ is plotted in the lower panel of figure 1 for our two fictitious stars; again the stationary star is the solid line and the receding star is the dotted line. When these terms are summed over $`j`$, we arrive at new determinations for the probability density of observing the star, indicated in the upper right panel of figure 1 by the horizontal lines. The original estimate of $`\varphi `$ was given by the middle dot and circle (at $`9\mathrm{kpc}`$), so it is evident that the “hill” and the “dip” in $`\varphi (r)`$ for the two stars are responsible for modifying the original estimate of $`\varphi `$ up or down in the correct direction.

The quantity $`\alpha `$ for the halo is usually taken to be between about 3 and 3.5 (Yanny & Newberg et al., 2000, and references therein). In the kinematic analysis below we assume $`\alpha =3.5`$, following the study of horizontal branch stars of Kinman, Suntzeff, & Kraft (1994). However, we find our results are fairly insensitive to $`\alpha `$; even $`\alpha =10`$ results only in small changes to our final results.

Another possible source of measurement error is in the radial velocities. As determined in Paper I, typical velocity errors may be as high as $`30\mathrm{km}\mathrm{s}^1`$. According to equation A3, the probability density of observing a star at a certain line-of-sight velocity $`v_{\mathrm{los}}`$ (with respect to the GC) follows a gaussian distribution, centered on $`v_{\mathrm{los}}=0\mathrm{km}\mathrm{s}^1`$, with a dispersion $`\widehat{\sigma }_{\mathrm{los}}`$ that depends on geometry and $`𝝈`$. Thus the probability density of observing the star at $`v_{\mathrm{los}}`$ considering a possible velocity error that follows a gaussian distribution with dispersion $`\mathrm{\Delta }v`$ is just

$`\varphi _\mathrm{e}(v_{\mathrm{los}})`$ $`=`$ $`{\displaystyle \frac{1}{\Gamma \sqrt{2\pi \mathrm{\Delta }v^2}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{v^2}{2\widehat{\sigma }_{\mathrm{los}}^2}}{\displaystyle \frac{(vv_{\mathrm{los}})^2}{2\mathrm{\Delta }v^2}}\right)𝑑v`$ (1)
$`=`$ $`{\displaystyle \frac{1}{\Gamma \sqrt{1+(\mathrm{\Delta }v/\widehat{\sigma }_{\mathrm{los}})^2}}}\mathrm{exp}\left({\displaystyle \frac{v_{\mathrm{los}}^2}{2(\widehat{\sigma }_{\mathrm{los}}^2+\mathrm{\Delta }v^2)}}\right);`$
note this reduces to the original equation A3 as $`\mathrm{\Delta }v0`$.

We have discussed the various modifications available to the maximum likelihood technique. In summary, they are the following:

* Subsample. Although we first determine the velocity ellipsoid using all stars, we can extract subsamples, specifically in certain magnitude bins, to investigate whether the velocity ellipsoid is a function of position in space.
* Spherical or cylindrical coordinates. The choice affects the values of $`𝐱`$ in the maximum likelihood procedure but is otherwise straightfoward.
* Magnitude errors. We either ignore magnitude errors or else use 11-point Gauss-Hermite integration with $`\mathrm{\Delta }m=0.5`$ and $`\alpha =3.5`$.
* Velocity errors. We either ignore velocity errors or else use equation 1 with $`\mathrm{\Delta }v=30\mathrm{km}\mathrm{s}^1`$.
* Number of fitting parameters. Several combinations of fitting parameters are used. We fit the three $`𝝈`$ parameters keeping $`(𝐯_{},v_{\mathrm{halo}})`$ at canonical values; we fit $`(v_{\mathrm{halo}},𝝈)`$ keeping $`𝐯_{}`$ at its canonical values; and we also fit all seven parameters simultaneously.

## 3 Results

Table 2 contains the determinations of $`(𝐯_{},v_{\mathrm{halo}},𝝈)`$ for several combinations of model assumptions discussed above; the first four columns describe these model assumptions. The first column states what criteria were used to select subsets of BHB stars from the total set of $`1170`$ BHB stars, and the number of stars in each subset is listed in the second column. The third column indicates whether a cylindrical or spherical geometry was assumed. The fourth column indicates how errors were handled; “0” indicates that errors were not considered. As described in §2, “m” indicates that 11-point Gauss-Hermite integration was used with $`\mathrm{\Delta }m=0.5`$, and “v” indicates that velocity errors were considered with $`\mathrm{\Delta }v=30`$. The following columns show the parameter determinations. All determinations in table 2 are for the three-parameter fit with $`(𝐯_{},v_{\mathrm{halo}})`$ fixed at canonical values (Dehnen & Binney, 1998) as shown, except for the second and third lines which fit four and seven parameters, respectively. The quoted errors consider correlations among the parameters, and the errors in parentheses do not consider the correlations. The former error is defined as $`\sqrt{W_i/2}`$ and the latter as $`1/\sqrt{2M_i}`$ where $`M_i`$ represents the diagonal elements of the second-derivative matrix of the log-likelihood function, and $`W_i`$ represents the diagonal elements of the inverse of that matrix.

The first row in Table 2 shows the results for the velocity ellipsoid for the complete sample ($`N_{\mathrm{stars}}=1170`$) in spherical coordinates, with no magnitude or velocity errors considered. The Sun’s velocity was fixed at the canonical value of $`𝐯_{}=𝐯_{\mathrm{LSR}}+(U,V,W)=(0,220,0)+(10,5,7)\mathrm{km}\mathrm{s}^1`$ (Dehnen & Binney, 1998). Under these assumptions all three components of the velocity ellipsoid are about $`100\mathrm{km}\mathrm{s}^1`$. The second row in Table 2 shows the results of simultaneously fitting $`v_{\mathrm{halo}}`$ and $`𝝈`$ to the same sample and otherwise same set of assumptions discussed above. The velocity ellipsoid has essentially the same determined values, and the rotation velocity of the halo is consistent with zero. The third row in Table 2 is the seven parameter fit to the sample. Again, the velocity ellipsoid is essentially the same, and the rotation velocity of the halo is seen to be marginally consistent with zero. Reassuringly, the velocity of the Sun with respect to the galactic halo in the $`y`$ direction (or equivalently, the Sun’s velocity in an inertial frame assuming the halo is actually non-rotating) is close to the expected value: $`v_{,y}v_{\mathrm{halo}}=222.1\pm 7.7\mathrm{km}\mathrm{s}^1`$. The error follows from rotating the eigenvectors by $`45\mathrm{deg}`$ and is smaller than either error for $`v_{,y}`$ or $`v_{\mathrm{halo}}`$ because the two quantities are strongly covariant. The $`x`$ and $`z`$ directions of $`𝐯_{}`$ are new measurements of the solar motion $`(U,W)`$. Our value of $`U`$ is consistent with that of Dehnen & Binney (1998); however our value of $`W`$ differs by several $`\sigma `$.

Figure 2 shows confidence regions for the components of the velocity ellipsoid for this case. Each panel in figure 2 is a two-dimensional parameter space where the other five parameters are held at their determined values and assumed to have symmetric errors, so the confidence regions are only approximatations to those one would obtain by integrating over the other five dimensions. Note that $`\sigma _r`$ is very well constrained. This is because the typical distance to the halo stars in our sample is large compared to the distance from the Sun to the galactic center, so measurements of $`v_{\mathrm{los}}`$ from Earth very nearly align with direct measurements of the galactocentric radial component of velocity. Also, each component of the velocity ellipsoid is anticorrelated with the other two; this demonstrates an intuitive “conservation law” of velocities. In the lower left panel, it is apparent that $`v_{\mathrm{halo}}`$ and $`v_{,y}`$ are correlated, so that the quantity $`v_{,y}v_{\mathrm{halo}}`$ is more tightly constrained than either quantity alone. This of course follows from the reasonable assumption that the halo rotation at the location of the Sun is in the $`y`$ direction.

The result that $`\sigma _r100\mathrm{km}\mathrm{s}^1`$ contrasts with most estimates of the velocity ellipsoid based on stellar samples in the solar neighborhood (Woolley, 1978; Pier, 1984; Carney & Latham, 1986; Norris, 1986; Layden et al., 1996). We now investigate some other sets of assumptions that may influence the results.

The fourth line in Table 2 introduces magnitude errors into the set of assumptions. These errors, at least in our implementation of a gaussian with width $`\mathrm{\Delta }m=0.5`$, apparently have little effect on the final results. Note that for a given kinematic model, distance affects the distribution of $`v_{\mathrm{los}}`$ only through the orientation of the velocity ellipsoid with respect to the line of sight. If the dispersions were isotropic and the halo nonrotating, which seem to be nearly the case, then $`v_{\mathrm{los}}`$ would be independent of distance. This may explain the insensitivity to magnitude errors. The fifth line in Table 2 incorporates both magnitude and velocity errors. Again, the final results are seen to remain relatively constant; however the values of the velocity ellipsoid components decrease modestly. In light of equation A3, the probability of observing a star at a certain velocity $`v_{\mathrm{los}}`$ is always higher for smaller values of $`|v_{\mathrm{los}}|`$; thus velocity errors tend to “favor” smaller velocities, which propagate to the smaller values of the velocity ellipsoid. However, the effect is small, of order $`\sqrt{\sigma ^2+\mathrm{\Delta }v^2}\sigma 5\mathrm{km}\mathrm{s}^1`$.

The next three lines in Table 2 show three-parameter determinations for the complete sample of BHB stars assuming the velocity ellipsoid is oriented along cylindrical coordinates. As for the spherical coordinates case, all three components of the velocity ellipsoid are still about $`100\mathrm{km}\mathrm{s}^1`$.

Determinations for cylindrical scenarios incorporating magnitude and velocity errors are shown in the tables, and as before, do not appear to affect the final results appreciably. Since velocity errors behave in a predictable way and have little effect on our final results, and magnitude errors have no noticable effect, we do not consider these errors further. It should be remembered that our error estimate was based on a set of highly simplified, and perhaps dubious, assumptions that do not take into account systematic errors, for instance.

In the following rows in Table 2 the sample of stars is modified. There are $`733`$ stars with $`g<18`$, so these stars extend from $`7\mathrm{kpc}`$ to $`30\mathrm{kpc}`$. The sample with $`g<16`$ contains the $`227`$ brightest stars (distances out to $`11.5\mathrm{kpc}`$). The velocity ellipsoid is seen to remain remarkably constant down to distances from the Sun only marginally outside the coverage of previous surveys (e.g., $`V13`$ for the sample of Layden et al. (1996)). The most intriguing result is that the radial component of the velocity ellipsoid $`\sigma _r`$ (or $`\sigma _R`$ in cylindrical coordinates) just outside the solar neighborhood is significantly smaller than the canonical value of $`150\mathrm{km}\mathrm{s}^1`$ in the solar neighborhood (Table 1). This incompatibility will be further discussed in §5.

The sample with $`g>18`$ contains the $`437`$ faintest stars. Since the stars are all very far away, the non-radial components of the velocity ellipsoid are very poorly constrained. Note that the quoted errors may also be ill-determined because of departures from gaussianity for small eigenvalues. However, the essential result remains the same: the radial component of the velocity ellipsoid is $`100\mathrm{km}\mathrm{s}^1`$.

Finally, to address concerns that the determination of the velocity ellipsoid may be contaminated by halo substructure, the last two lines of Table 2 show the results for all stars excluding the Sagittarius stream. This exclusion is conservatively accomplished by simply rejecting all BHB stars within 10 degrees of the great circle with pole $`(\mathrm{},b)=(94,11)`$, which approximates the Sagittarius stream (Ivezić et al., 2003, Paper I). Since the Sagittarius stream is the most conspicuous substructure in our sample of BHB stars (Paper I), the stability of the velocity ellipsoid solution indicates that halo substructure does not appreciably influence the results.

## 4 Monte Carlo simulations

An independent way of verifying the accuracy of the maximum likelihood technique is through Monte Carlo simulations. In case the geometry of the SDSS footprint (see Paper I) affects the analysis, we assemble $`1170`$ artificial stars with the same angular coordinates and measured $`g`$ magnitudes as the complete BHB sample, but with newly constructed velocities. For simplicity absolute magnitudes are assumed constant at $`0.7`$. Only velocity ellipsoids aligned with spherical coordinates will be considered here. These Monte Carlo simulations perform the three-parameter fit to the velocity ellipsoid with $`𝐯_{}=(10,225,7)\mathrm{km}\mathrm{s}^1`$ and $`v_{\mathrm{halo}}=0\mathrm{km}\mathrm{s}^1`$. Note, however, that the values of $`𝐯_{}`$ and $`v_{\mathrm{halo}}`$ are irrelevant here since they are first added to the stars’ velocities in construction, and then subtracted as known in the parameter determination.

We run 100 Monte Carlo simulations on each of four setups. The first setup does not consider magnitude or velocity errors; the second setup allows for $`\mathrm{\Delta }m=0.5`$; the third assumes $`\mathrm{\Delta }v=30\mathrm{km}\mathrm{s}^1`$. The fourth setup considers the magnitude and velocity errors as well as allowing for contamination.

The constructed set of artificial stars has known space positions, so it is straightforward to apply an assumed velocity ellipsoid to obtain a simulated three-dimensional space velocity for each star. Here we assume $`𝝈=(150,100,80)\mathrm{km}\mathrm{s}^1`$ to strain the maximum likelihood analysis. The three-dimensional velocities are then projected onto the line of sight to each star to obtain $`v_{\mathrm{los}}`$, and then the projections of $`𝐯_{}`$ and $`v_{\mathrm{halo}}`$ are subtracted to obtain $`v_r`$. This comprises the simulated data set for the first setup. To incorporate magnitude errors, the radial distances to the stars are shuffled, i.e. the $`g`$ magnitudes are shifted by a gaussian error distribution of width $`\mathrm{\Delta }m`$, $`v_r`$ is calculated as before, and then the $`g`$ magnitudes are shifted back to their “observed” positions. Velocity errors are incorporated by adding an error to each observed $`v_r`$. To incorporate contamination, a fraction of stars are randomly flagged as disk stars. These stars are then assumed to be 2.5 times closer (2 mag) to the Sun than they would be if they were BHB stars. Their velocities are taken to be $`220\mathrm{km}\mathrm{s}^1`$ in the $`\varphi `$ direction plus three random components from a *disk* velocity dispersion, which is assumed isotropic with components $`(20,20,20)\mathrm{km}\mathrm{s}^1`$. Dehnen & Binney (1998) find a non-isotropic disk velocity dispersion with components of the same order (depending on color), but our intent is to contaminate the simulated halo star population in a simple and controlled way.

Figure 3 shows the results of these Monte Carlo simulations with the four setups colored black (no errors), blue (magnitude errors), red (velocity errors), and green (magnitude and velocity errors with 10% contamination). As shown, the scatter in the results of the Monte Carlo simulations is comparable to the confidence regions shown in figure 2, so our error analysis is verified to be reasonable. As in the analysis of the real data, the radial component of the velocity ellipsoid $`\sigma _r`$ is very well constrained; one can also see anticorrelation between any two components. The black dots and blue dots occupy essentially the same region of parameter space; this verifies, independently of the method used in §3, that magnitude errors of the order $`\mathrm{\Delta }m`$ do not affect our results significantly.

As shown in Table 3, the mean determinations of $`𝝈`$ over 100 Monte Carlo simulations recover the input values for the “no errors” case and the “magnitude errors only” case. It is also clear that velocity errors of scale $`\mathrm{\Delta }v=30\mathrm{km}\mathrm{s}^1`$ push the measured values higher by $`5\mathrm{km}\mathrm{s}^1`$ because the error adds in quadrature. This result is consistent with the independent determination in §3 which stated that with $`\mathrm{\Delta }v=30\mathrm{km}\mathrm{s}^1`$ the “true” values of the velocity ellipsoid (i.e., those measured using the analysis that allows for velocity errors) were $`5\mathrm{km}\mathrm{s}^1`$ lower than the “measured” values (not considering velocity errors). Finally, it is interesting that for the case which includes contamination by disk stars, the determination of $`\sigma _\varphi `$ increases while those of $`\sigma _r`$ and $`\sigma _\theta `$ both decrease. This is because the disk stars are constructed to travel primarily in the $`\varphi `$ direction.

However, the main conclusion to be gained from these Monte Carlo simulations is that the radial component of the velocity ellipsoid never strays too far from its correct value of $`150\mathrm{km}\mathrm{s}^1`$, even allowing for errors such as contamination. Thus the maximum likelihood analysis is verified to be reliable and stable.

## 5 Summary and discussion

Applying maximum-likelihood methods to the line-of-sight velocities and photometric distances of BHB stars, we have modeled the kinematics of the relatively distant galactic halo. (The median distance of the sample is about $`25\mathrm{kpc}`$.) No significant anisotropy of the velocity ellipsoid is detected. In a model aligned with spherical coordinates, for example, $`𝝈=(101.4,97.7,107.4)\mathrm{km}\mathrm{s}^1`$. However, the errors in $`\sigma _\theta `$ and $`\sigma _\varphi `$ are much larger, $`1520\mathrm{km}\mathrm{s}^1`$, than that of $`\sigma _r`$, because of the great distances of most stars: that is, the Galactocentric radius $`r`$ of most stars is more than twice that of the Sun, so that the direction of the line of sight is not far from radial despite the broad angular range of our sample. This result contrasts starkly and surprisingly with the radial anisotropy of halo velocities in the solar neighborhood ($`𝝈(150,100,100)\mathrm{km}\mathrm{s}^1`$; see Table 1). Our results are consistent, however, with the only other comparable study of distant halo field stars: Sommer-Larsen et al. (1997).

We are able to fit separately for the rotation speed of the halo ($`v_{\mathrm{halo}}`$) and the orbital velocity of the Sun, although the two are strongly covariant. The halo rotation is marginally consistent with zero at the one-sigma level, which is $`20\mathrm{km}\mathrm{s}^1`$. The difference $`v_{,y}v_{\mathrm{halo}}=222.1\pm 7.7\mathrm{km}\mathrm{s}^1`$ is more tightly constrained than either quantity separately. This is consistent with a non-rotating halo if one accepts the IAU recommendation for the local standard of rest, $`𝐯_{\mathrm{LSR}}=220\mathrm{km}\mathrm{s}^1`$ (Kerr & Lynden-Bell, 1986) and Dehnen & Binney (1998)’s solar motion, $`𝐯_{}𝐯_{\mathrm{LSR}}=(10,5,7)\mathrm{km}\mathrm{s}^1`$, although some authors have advocated a substantially smaller (Olling & Merrifield, 1998) or larger (Uemura et al., 2000) value of the azimuthal component of $`𝐯_{\mathrm{LSR}}`$. (If the solar rotation *is* substantially less than $`225\mathrm{km}\mathrm{s}^1`$, then according to our results, the halo should counterrotate!) Our result is also statistically consistent with most of the studies cited in Table 1, although the latter seem to prefer a slight positive halo rotation. Note that what those studies actually measure is the mean velocity of local halo stars relative to the Sun; the values for $`v_{halo}`$ in the fourth column of that table have been adjusted to $`V_{}=225\mathrm{km}\mathrm{s}^1`$ so that they can be compared directly to our result. Presumably, $`𝐯_{\mathrm{LSR}}`$ will soon be determined with exquisite accuracy by a variant of the classical visual-binary method applied to stars orbiting Sgr $`\mathrm{A}^{}`$, but accurate values are not in the public domain at the time of writing (Ghez et al., 2003; Salim & Gould, 1999).

Of course, the measurement of $`v_{\mathrm{halo}}`$ depends upon the assumption that the halo rotation curve is flat; line-of-sight velocities would be unaffected by adding a solid-body component to the angular velocity of the halo, although proper motions would be (see below).

These are the main results of our study. It is likely that more science could be extracted from this sample, and completion of the SDSS will allow it to be enlarged. For the present, we simply highlight certain interesting features of the data that have not yet been thoroughly explored.

The distribution of $`v_{\mathrm{los}}`$ is plotted in figure 4 and beautifully follows a gaussian of width $`\sigma _{\mathrm{los}}=101.6\pm 3.0\mathrm{km}\mathrm{s}^1`$. The uncertainty quoted here is one-sigma statistical, to which should be added a similar systematic uncertainty pending better quantification of our measurement errors. An estimate of this systematic error is easily determined as $`\sigma _{\mathrm{los}}\sqrt{\sigma _{\mathrm{los}}^2\mathrm{\Delta }v^2}3.4\mathrm{km}\mathrm{s}^1`$ ($`\mathrm{\Delta }v=26\mathrm{km}\mathrm{s}^1`$ from Paper I). The line-of-sight velocities (figure 4) are about as gaussianly distributed as the number of stars would allow. This is consistent with (although not required by) the isotropy of the velocity ellipsoid that we find, since a strongly anisotropic ellipsoid, even if gaussian along each individual axis, would not produce an exact gaussian in an average over many lines of sight. Gaussianity could result from large observational errors, regardless of the intrinsic velocity distribution, if the errors were gaussian, but our errors seem not to be large enough to dominate the observed distribution (Paper I). By contrast, the distribution of transverse velocities of halo stars in the solar neighborhood is distinctly nongaussian, with a negative kurtosis (Popowski & Gould, 1998). This is another indication that the local halo differs from the more distant halo probed by our BHB sample.

Another aspect of the data that could be thoroughly explored is the smoothness of the BHB phase-space distribution. The notion of a stellar velocity ellipsoid implicitly presumes that phase space is smoothly populated. If instead the distribution is very lumpy and the number of important lumps or streams is small, then statistical characterizations of *individual* halos, in particular that of the Galaxy, may not be meaningful. In cosmology, however, statistical descriptions are unavoidable. Current theory holds that galactic halos have been built up by accretion of lesser systems (Searle & Zinn, 1978; Bullock, Kravtsov, & Weinberg, 2001). Accreted gas joins galactic disks, and the details of its history are effaced by dissipation. Stars that form in these systems before they are accreted (and presumably also their dark matter), however, behave collisionlessly and preserve a partial record of their origins. The structure of the stellar and dark-matter halos should be related. Helmi, White, & Springel (2003) conclude from a detailed analysis and scaling of cosmological simulations that the number of streams composing a halo similar to that of the Galaxy should be large, amounting to almost $`10^4`$ within $`25\mathrm{kpc}`$, and that most streams joined the Galaxy early on, allowing more time for phase mixing. The number of our BHB stars is much smaller than this, so that one might expect to find a smooth distribution. And yet we see structure in our sample that is consistent with the tidal stream of the Sagittarius Dwarf (Paper I), as previously delineated in SDSS data (Ivezić et al., 2000; Yanny & Newberg et al., 2000; Ibata et al., 2001a; Newberg & Yanny et al., 2002) and in a survey of carbon stars (Ibata et al., 2001b). Most of our results for the halo velocity ellipsoid are scarcely affected by excluding that part of the sample that lies close to the Sagittarius stream on the sky. In short, it is not yet clear whether the kinematics of the Galactic halo admit a statistical decription. Larger samples of stars will be necessary to clarify this.

In Paper I, we used proper motions only as a check on the purity of our sample. Calibration against the QSOs shows that the proper motions of individual BHB stars in our sample are smaller than measurement errors by factors $`35`$, but statistical detection of the kinematics of the population as a whole should certainly be possible. Among other benefits, this would allow a direct constraint on the solid-body component of halo rotation and a further test of the isotropy of the velocity dispersions. Systematic errors in the proper motion data present many complexities, however, so we defer their analysis to a later paper. The quality of the USNO-B data should be appreciably better than that of the USNO-A, and one should probably use the former to determine the proper motions of SDSS objects.

Our BHB radial velocity data could be used to constrain the mass of the outer Galaxy. We have not yet attempted a serious dynamical analysis, partly in order to do full justice to the kinematics, but also because mass measurements rely more heavily on aspects of the data that are less certain: namely, the photometric distances and the completeness as a function of magnitude. However, since no kinematic information was used in the selection of the sample, we may assume that the velocities fairly represent the whole population of BHB stars, whose spatial density is often described as a powerlaw in Galactocentric radius, $`nr^\alpha `$, with exponent $`\alpha 3.5`$ (Kinman, Suntzeff, & Kraft, 1994). In that case, to the extent that the dispersions are isotropic and constant with radius and rotation is negligible, and that the potential is spherical, Jeans’ equations imply

$`v_{\mathrm{circ}}`$ $``$ $`187\left({\displaystyle \frac{\alpha }{3.5}}\right)^{1/2}\left({\displaystyle \frac{\sigma }{100\mathrm{km}\mathrm{s}^1}}\right),`$ (2)
$`M(r)`$ $``$ $`2.0\times 10^{11}\left({\displaystyle \frac{v_{\mathrm{circ}}}{187\mathrm{km}\mathrm{s}^1}}\right)^2\left({\displaystyle \frac{r}{25\mathrm{kpc}}}\right)M_{}.`$ (3)
The value of $`v_{\mathrm{circ}}`$ so obtained, with $`\sigma 100\mathrm{km}\mathrm{s}^1`$, is substantially lower than the circular velocity in the vicinity of the Sun. Faced with similar evidence drawn from their own independent BHB sample, Sommer-Larsen et al. (1997) concluded that the tangential velocity dispersion of the halo probably rises beyond $`r20\mathrm{kpc}`$ so as to exceed the radial dispersion and be consistent with a constant circular velocity. At such distances, however, the tangential dispersions are very poorly constrained by line-of-sight velocities without *a priori* assumptions about the potential (see Table 2). Rather than make such assumptions, one may prefer to await further direct evidence.

We thank A. Gould, B. Paczyński, R. Lupton, and J. Hennawi for helpful discussions, and M. Strauss for a close reading of a draft of this paper. GRK is grateful for generous research support from Princeton University and from NASA via grants NAG-6734 and NAG5-8083. Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, the University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington.

## Appendix A The maximum likelihood technique

In this section we outline the maximum likelihood procedure used to determine the velocity ellipsoid of the Galactic halo.

Halo stars are assumed to have a gaussian velocity distribution whose principal axes align with spherical or cylindrical coordinates. In this appendix we assume the velocity of the Sun relative to the galactic center is known and subtracted from stars’ radial velocities, and that every star’s position in space is known exactly; see §2 for a discussion of magnitude and velocity errors. We define $`𝐱`$ as the normalized heliocentric line-of-sight vector to the star, so that the line-of-sight velocity $`v_{\mathrm{los}}=𝐱𝐯`$ where $`𝐯`$ is the true space velocity of the star. The three components of each of the vectors $`𝐱`$, $`𝐯`$, and $`𝝈`$ should be thought of as projections along the principal axes of the velocity ellipsoid at the specific point in space. The probability $`\varphi dv_{\mathrm{los}}`$ of observing a given star with a radial velocity between $`v_{\mathrm{los}}`$ and $`dv_{\mathrm{los}}`$ is expressed as an integral over the possible combination of true space velocities that would result in an observation between $`v_{\mathrm{los}}`$ and $`dv_{\mathrm{los}}`$:

$$\varphi (v_{\mathrm{los}})dv_{\mathrm{los}}=_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\varphi _1(v_1)\varphi _2(v_2)\varphi _3(v_3)𝑑v_1𝑑v_2𝑑v_3$$
(A1)
where $`v_1`$, $`v_2`$, and $`v_3`$ are the components of the star’s velocity in the directions of the principal axes. Here $`\varphi _i(v_i)dv=(2\pi \sigma _i^2)^{1/2}\mathrm{exp}(v_i^2/2\sigma _i^2)dv`$ where $`\sigma _i`$ is the component of the velocity ellipsoid in the $`i`$ direction. Geometry dictates a relation between $`dv_{\mathrm{los}}`$ and the other infinitesimal quanitites: $`dv_{\mathrm{los}}=|x_i|dv_i`$. Also, the vector $`𝐯`$ is constrained by $`v_{\mathrm{los}}`$ by the relation $`v_{\mathrm{los}}=𝐱𝐯`$ so finally our integral becomes

$$\varphi (v_{\mathrm{los}})dv_{\mathrm{los}}=\frac{dv_{\mathrm{los}}}{|x_1|}_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑑v_2𝑑v_3\varphi _1\left(\frac{1}{x_1}(v_{\mathrm{los}}x_2v_2x_3v_3)\right)\varphi _2(v_2)\varphi _3(v_3)$$
(A2)
Here the integral integrates over the $`2`$ and $`3`$ directions, which is mathematically correct if $`x_10`$, but to avoid small-number numerical problems, in practice the code shuffles the indices as appropriate to ensure that the largest component of $`𝐱`$ is placed in the denominator. The value of the integral is

$`\varphi (v_{\mathrm{los}})`$ $`=`$ $`{\displaystyle \frac{1}{\Gamma }}\mathrm{exp}\left({\displaystyle \frac{v_{\mathrm{los}}^2}{2\widehat{\sigma }_{\mathrm{los}}^2}}\right)`$ (A3)
$`{\displaystyle \frac{1}{2\widehat{\sigma }_{\mathrm{los}}^2}}=C^{}B^2/4A^{}\Gamma `$ $`=`$ $`|x_1|\sigma _1\sigma _2\sigma _3\sqrt{8\pi AA^{}}`$
$`A^{}={\displaystyle \frac{x_3^2}{2\sigma _1^2x_1^2}}(D1)+{\displaystyle \frac{1}{2\sigma _3^2}}B^{}`$ $`=`$ $`{\displaystyle \frac{x_3}{\sigma _1^2x_1^2}}(D1)C^{}={\displaystyle \frac{1}{2\sigma _1^2x_1^2}}(D1)`$
$`A={\displaystyle \frac{1}{2\sigma _2^2}}+{\displaystyle \frac{x_2^2}{2\sigma _1^2x_1^2}}D`$ $`=`$ $`{\displaystyle \frac{\sigma _2^2x_2^2}{\sigma _1^2x_1^2+\sigma _2^2x_2^2}}.`$
Note that at a given position, stars are most likely to be found at zero velocity (with respect to the galactic center) but have a characteristic velocity dispersion $`\widehat{\sigma }_{\mathrm{los}}`$, which is a function of the geometry of the stars with respect to the Sun and of $`𝝈`$. Since the probability of observing a set of stars with known space locations and observed line-of-sight velocities is the product of the individual probabilities, the log-likelihood is

$$(𝝈)=\frac{N}{2}\mathrm{ln}\left(8\pi \sigma _1^2\sigma _2^2\sigma _3^2\right)+\left[\left(\frac{B^2}{4A^{}}C^{}\right)v_{\mathrm{los}}^2\frac{1}{2}\mathrm{ln}(x_1^2AA^{})\right]$$
(A4)
where the summation is over all observed stars. Assuming a uniform prior (Lupton, 1993), we find the maximum of $`(\sigma _1,\sigma _2,\sigma _3)`$ to determine the most likely values of the velocity ellipsoid $`𝝈`$.

If quantities such as the velocity of the Sun with respect to the Galactic center or the bulk rotation velocity of the halo are free parameters, they can be determined from this maximum-likelihood estimate. These quantities affect $`v_{\mathrm{los}}`$ for every star, so it is straightforward to extrapolate equation A4 to $`=(𝝈,𝐯_{},v_{\mathrm{halo}})`$.


format: text

text:
# Some exact results for the trapping of subdiffusive particles in one dimension

## 1 Introduction

The trapping of Brownian particles by static traps randomly distributed over either a Euclidean or a disordered substrate is a fundamental problem of non-equilibrium statistical mechanics and chemistry with a very wide range of applications . This is also one of the oldest problems in random walk theory which was essentially formulated by Smoluchowski at the beginning of past century in his theory of coagulation of colloidal particles . This model has proven useful in research areas such as the trapping of mobile defects in crystals with point sinks , the kinetics of luminescent organic materials , the kinetics of photosynthetic light energy to oxygen conversion , anchoring of polymers by chemically active sites , atomic diffusion in glasslike materials , and others . Some generalizations of the basic model have also been considered recently. The case of a gated (ungated) random walker trapped by a distribution of ungated (gated) fixed traps has been proposed in connection with the kinetics of reactions between complex molecules whose active groups are screened as the molecules diffuse . The trapping problem with many random walkers has also been studied in one dimension  and the corresponding trapping statistics have been found on Euclidean and fractal lattices . The multiparticle predator-prey problems in which a single fixed trap, the “lamb”, is captured by one of a set of $`N`$ random walkers or “lions” initially placed at a given distance from the prey was also discussed by Krapivsky and Redner . Another quite interesting and difficult variation of the standard trapping problem is that where the traps diffuse too .

Trapping reactions between molecules embedded in biological samples and disordered materials are usually handicapped by the porous and statistical fractal structure of these media . In some cases this gives rise to subdiffusion of the particles, i.e., the mean square displacement $`r^2(t)`$ of the particles from the original starting site is no longer linear on time, but verifies a generalized Fick’s second law:

$$r^2(t)\frac{2K_\gamma }{\mathrm{\Gamma }(1+\gamma )}t^\gamma $$
(1)
where $`\gamma `$ (with $`0<\gamma <1`$) is the (anomalous) diffusion exponent and $`K_\gamma `$ is the diffusion coefficient. Of course, there are many other instances in which subdiffusion processes appear . A useful approach for understanding subdiffusion processes is by means of the continuous time random walk model (CTRW) in which the random walker performs jumps with a waiting time distribution with a broad long-time tail: $`\psi (t)t^{(1+\gamma )}`$ for large $`t`$ . The long-time tail of this waiting time distribution incorporates, in a statistical sense, the effect of the bottlenecks and the dead-ends in the diffusion of the random walker embedded in the disordered structure, and the model is compatible with Eq. (1). For subdiffusive random walkers the continuum description given by the ordinary diffusion equation is replaced by the fractional diffusion equation

$$\frac{}{t}W(x,t)=K_\gamma {}_{0}{}^{}D_{t}^{1\gamma }\frac{^2}{x^2}W(x,t),$$
(2)
where $`{}_{0}{}^{}D_{t}^{1\gamma }`$ is the Riemann-Liouville fractional derivative of order $`1\gamma `$ ,

$${}_{0}{}^{}D_{t}^{1\gamma }W(x,t)=\frac{1}{\mathrm{\Gamma }(\gamma )}\frac{}{t}_0^t𝑑\tau \frac{W(x,\tau )}{(t\tau )^{1\gamma }},$$
(3)
and $`W(x,t)`$ is the probability density that the particle that started at $`0`$ at time $`0`$ is at $`x`$ at time $`t`$.

Our objective in this paper is to find analytical expressions for the survival probability $`\mathrm{\Phi }(t)`$ defined as the probability that no trap site has been reached by the *subdiffusive* random walker by time $`t`$. Two variations of the classical trapping problem (sometimes called Rosenstock’s trapping problem) are considered: (i) the “one-sided” trapping problem  in which only one half-line of a one-dimensional lattice is filled with a random distribution of static traps (this process could mimic the excitation or production of defects on one side of a fiber by irradiation, the other side being shielded); and (ii) the “two-sided” trapping problem corresponding to the trapping of a single subdiffusive random walker placed initially at $`x=0`$ between two half-lines ($`x<0`$ and $`x>0`$) randomly filled with static traps. We find exact analytical solutions for both problems and compare them with simulations. As precedents of these results we whould cite the work of Blumen, Klafter and Zumofen (see  and references therein) about the trapping of particles on Euclidean and fractal substrates for random walkers with waiting time distributions with broad long-time tails: $`\psi (t)t^{1\gamma }`$, $`0<\gamma <1`$.

The paper is organized as follows. In Sec. 2 the exact solution of the one-sided trapping problem is obtained and compared with simulation results. In Sec. 3 we build the extended Rosenstock’s approximation for the survival probability of the particle in the one-sided trapping problem by calculating an exact expression for the moments of the one-sided span (the number of distinct sites visited by the random walker in a given direction). In Sec. 4 we obtain the exact survival probability of the subdiffusive random walker for the two-sided case in integral form. This integral is evaluated in the asymptotic long-time and short-time limits in order to get simple closed expressions. The paper ends with some conclusions and remarks in Sec. 5.

## 2 The one-sided trapping problem

In the one-sided trapping model quenched traps are randomly distributed on, say, the right-hand side of a one-dimensional lattice ($`x>0`$) with concentration $`c`$. The random walker is placed initially upon site $`x=0`$ and performs jumps to its nearest neighbour sites with a waiting time distribution $`\psi (t)`$ until it reaches a trap site where it is absorbed. The survival probability is given by :

$$\mathrm{\Phi }(t)=\underset{r=1}{\overset{t}{}}e^{\lambda r}P(t|r),$$
(4)
where $`\lambda =\mathrm{ln}(1c)`$ and $`P(t|r)`$ is the probability that the span of the random walker in the positive direction (the largest distance reached by the random walker for $`x>0`$) is equal to $`r`$ after $`t`$ time steps. Now let $`\mathrm{\Gamma }(t|r)`$ be the probability that the site $`x=r`$ has not been visited by the random walker by time $`t`$ (the so-called fixed-trap survival probability). Because

$$P(t|r)=\frac{d\mathrm{\Gamma }(t|r)}{dr},$$
(5)
Eq. (4) takes the following form in the continuous limit:

$$\mathrm{\Phi }(t)=_0^{\mathrm{}}e^{\lambda r}\frac{d\mathrm{\Gamma }(t|r)}{dr}𝑑r,$$
(6)
or, integrating by parts,

$$\mathrm{\Phi }(t)=\lambda _0^{\mathrm{}}e^{\lambda r}\mathrm{\Gamma }(t|r)𝑑r,$$
(7)
because $`\mathrm{\Gamma }(t|0)=0`$, $`\mathrm{\Gamma }(t|\mathrm{})=1`$ and $`e^{\lambda r}0`$ as $`r\mathrm{}`$. For subdiffusive particles, the function $`\mathrm{\Gamma }(t|r)`$ can be written in terms of Fox’s $`H`$ function :

$$\mathrm{\Gamma }(t|r)=1H_{11}^{10}\left[\frac{r}{\sqrt{K_\gamma }t^{\gamma /2}}|\begin{array}{c}(1,\gamma /2)\hfill \\ (0,1)\hfill \end{array}\right].$$
(8)
The time Laplace transform of $`\mathrm{\Gamma }(t|r)`$ is given by

$$\stackrel{~}{\mathrm{\Gamma }}(s|r)=\frac{1}{s}\left[1\text{exp}(r\sqrt{s^\gamma })\right].$$
(9)
Then the Laplace transform of $`\mathrm{\Phi }(t)`$ is readily calculated from Eqs. (7) and (9):

$$\stackrel{~}{\mathrm{\Phi }}(s)=\frac{1}{s+\lambda \sqrt{K_\gamma }s^{1\gamma /2}}$$
(10)
so that

$$\mathrm{\Phi }(t)=E_{\gamma /2}(\xi ),$$
(11)
where $`\xi \lambda \sqrt{K_\gamma t^\gamma }`$ and $`E_\alpha (z)`$ is the Mittag-Leffler function with parameter $`\alpha `$ . For $`\gamma =1`$, the Mitagg-Leffler function becomes :

$$\mathrm{\Phi }(t)=e^{\xi ^2}\text{erfc}\left(\xi \right),$$
(12)
with $`\xi =\lambda \sqrt{Dt}`$, ($`K_1D`$), and we recover the result for the survival probability of a normal diffusive random walker in the one-sided trapping problem . For very long times, the asymptotic expansion of the Mittag-Leffler function  allows us to write

$$\mathrm{\Phi }(t)=\underset{k=1}{\overset{n}{}}\frac{(1)^{k+1}}{\mathrm{\Gamma }(1k\gamma /2)}\xi ^k+𝒪(\xi ^{1n}).$$
(13)
Thus, an asymptotic time regime is reached for $`t1/(K_\gamma \lambda ^2)^{1/\gamma }`$ where the survival probability exhibits a power law decay

$$\mathrm{\Phi }(t)\frac{1}{\mathrm{\Gamma }(1\gamma /2)}\frac{1}{\lambda \sqrt{K_\gamma t^\gamma }}.$$
(14)
This is an algebraic fluctuation slowdown corresponding to the Donsker-Varadhan limit . We compare in Fig. 1 the exact survival probability for a subdiffusive random walker with a waiting time distribution with the form of the Pareto law

$$\psi (t)=\frac{\gamma }{(1+t)^{1+\gamma }}$$
(15)
and $`\gamma =1/2`$. Taking the distance between nearest neighbour sites (jump length) as 1, the subdiffusion constant is $`K_{1/2}=1/\left[2\mathrm{\Gamma }(1/2)\right]=1/\sqrt{4\pi }`$. The concentration of traps used was $`c=0.01`$.

## 3 The Rosenstock’s approximation for the one-sided trapping problem

Let $`S(t)`$ be the number of distinct sites on the positive half-line visited up to time $`t`$ by a random walker who started at $`x=0`$ at time $`t=0`$. Then, the survival probability of the random walker is given by $`\mathrm{\Phi }(t)=(1c)^{S(t)}=e^{\lambda S(t)}`$, the average being performed over all realizations of the random walker’s exploration of the lattice from time $`0`$ until time $`t`$. By means of the cumulant expansion technique  the following equivalent form is derived:

$$\mathrm{\Phi }(t)=\mathrm{exp}\left[\underset{n=1}{\overset{\mathrm{}}{}}(1)^n\frac{\kappa _n\lambda ^n}{n!}\right],$$
(16)
where $`\kappa _n`$, $`n=1,2,\mathrm{}`$ denote the cumulants of $`S(t)`$: $`\kappa _1=S(t)`$, $`\kappa _2=S^2(t)S(t)^2\sigma ^2(t)`$, $`\mathrm{}`$. If we keep the first $`n+1`$ terms of the sum in equation (16) we arrive at the $`n`$th-order Rosenstock approximation . The error made by using this approximation is $`𝒪(\kappa _{n+2}\lambda ^{n+2})`$. Thus, the condition $`\kappa _{n+2}1/\lambda ^{n+2}`$ must be fulfilled for the $`n`$th-order Rosenstock approximation to be reasonable.

For the Rosenstock approximation to be useful for the subdiffusive case, it is necessary to know the cumulants or, equivalently, the moments of $`S(t)`$ for subdiffusive particles. The next step is thus to derive exact expressions for the moments $`S^m(t)`$, $`m=1,2,\mathrm{}`$ of the one-sided span of a subdiffusive random walker, i.e., the territory explored in a given direction. The probability density of the one-sided span, $`P(t|r)`$, was given in Eq. (5) in terms of the survival probability. Hence

$$S^m(t)=_0^{\mathrm{}}\frac{d\mathrm{\Gamma }(t|r)}{dr}r^m𝑑r=_0^{\mathrm{}}\frac{d}{dr}(1\mathrm{\Gamma }(t|r))r^m𝑑r$$
(17)
or, integrating by parts,

$$S^m(t)=m_0^{\mathrm{}}\left[1\mathrm{\Gamma }(t|r)\right]r^{m1}𝑑r,$$
(18)
where we have taken into account that $`\mathrm{\Gamma }(t|0)=0`$ and $`lim_r\mathrm{}r^m[1\mathrm{\Gamma }(t|r)]=0`$ as a consequence of the stretched exponential behaviour of $`1\mathrm{\Gamma }(t|r)`$ as $`r\mathrm{}`$ for subdiffusive particles . From Eq. (18) and the explicit expression of the survival probability in Eq. (8) we obtain

$$S^m(t)=ms_m\left(K_\gamma t^\gamma \right)^{m/2},$$
(19)
with

$$s_m=_0^{\mathrm{}}z^{m1}H_{11}^{10}\left[z|\begin{array}{c}(1,\frac{\gamma }{2})\hfill \\ (0,1)\hfill \end{array}\right]𝑑z.$$
(20)
In order to calculate the coefficients $`s_m`$, $`m=1,2,\mathrm{}`$ we observe that the Laplace transform of the Fox $`H`$ function appearing in the integrand in Eq. (20), $`\stackrel{~}{H}(u)`$, is a generating function of these coefficients:

$$\stackrel{~}{H}(u)=\underset{m=1}{\overset{\mathrm{}}{}}(1)^{m1}\frac{u^{m1}}{(m1)!}s_m.$$
(21)
Taking into account some properties of the Fox functions  we finally identify $`\stackrel{~}{H}(u)`$ with a two-parameter Mittag-Leffler function as follows:

$$\stackrel{~}{H}(u)=E_{\gamma /2,1+\gamma /2}(u),$$
(22)
where, in the last identity, we have used the relation

$$E_{\alpha \beta }(u)=H_{12}^{11}\left[u|\begin{array}{c}(0,1)\hfill \\ (0,1),(1\beta ,\alpha )\hfill \end{array}\right]$$
(23)
between the two-parameter Mittag-Leffler functions and the Fox functions . The function $`E_{\alpha \beta }(u)`$ admits the following series expansion :

$$E_{\alpha \beta }(u)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{z^k}{\mathrm{\Gamma }(\alpha k+\beta )},\alpha >0,\beta >0.$$
(24)
From Eqs. (21), (22) and (24) we find

$$s_m=\frac{(m1)!}{\mathrm{\Gamma }(1+m\gamma /2)},m=1,2,\mathrm{}$$
(25)
so that from Eq. (19) we finally obtain

$$S^m(t)=\frac{m!}{\mathrm{\Gamma }(1+m\gamma /2)}(K_\gamma t^\gamma )^{m/2}.$$
(26)
Obviously, for $`\gamma =1`$ we recover the exact result for the moments of the one-sided span of a diffusive random walker :

$$S^m(t)=\frac{\mathrm{\Gamma }[(m+1)/2]}{\sqrt{\pi }}(4Dt)^{m/2}.$$
(27)
From Eq. (26) one can calculate the cumulants $`\kappa _n`$, $`n=1`$,$`2`$,$`\mathrm{}`$ of the one-sided span distribution, and their direct substitution into Eq. (16) yields the general expression for the Rosenstock approximation for subdiffusive particles:

$$\mathrm{\Phi }(t)=\mathrm{exp}\left\{\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1)^n}{n!}a_n\xi ^n\right\}.$$
(28)
The first three coefficients $`a_n`$ are:

$$\begin{array}{ccc}\hfill a_1& =& \frac{1}{\mathrm{\Gamma }(1+\frac{\gamma }{2})}\hfill \\ \multicolumn{3}{c}{}\\ \hfill a_2& =& \frac{1}{[\mathrm{\Gamma }(1+\frac{\gamma }{2})]^2}+\frac{2}{\mathrm{\Gamma }(1+\gamma )}\hfill \\ \multicolumn{3}{c}{}\\ \hfill a_3& =& \frac{2}{[\mathrm{\Gamma }(1+\frac{\gamma }{2})]^3}\frac{6}{\mathrm{\Gamma }(1+\frac{\gamma }{2})\mathrm{\Gamma }(1+\gamma )}+\frac{6}{\mathrm{\Gamma }(1+\frac{3\gamma }{2})}.\hfill \end{array}$$
(29)
It is interesting to note that from Eqs. (11) and (28) we get the following identity

$$\mathrm{ln}E_{\gamma /2}(\xi )=\underset{n=1}{\overset{\mathrm{}}{}}\frac{a_n}{n!}\xi ^n,$$
(30)
which provides a series expansion of the logarithm of the Mittag-Leffler function.

## 4 The two-sided trapping problem

In this section we calculate the exact survival probability of a single subdiffusive random walker placed initially between two half-lines populated with a random distribution of traps with concentration $`c`$. We will use parallel arguments to the standard ones for diffusive random walkers . We must also mention that an alternative approach was proposed by Anlauf  to find the long-time behavior for the one-dimensional trapping problem with a diffusive particle.

As a starting point we evaluate the probability $`W(x,t|x_0,t=0)`$ that a random walker starting from $`x_0`$ at $`t=0`$ inside a box $`0xL`$ with absorbing boundaries is at $`x`$ at time $`t`$. The function $`W(x,t|x_0,t=0)`$ satisfies the fractional partial differential equation (the subdiffusion equation):

$$\frac{}{t}W(x,t|x_0,0)=K_\gamma {}_{0}{}^{}D_{t}^{1\gamma }\frac{^2}{x^2}W(x,t|x_0,0),$$
(31)
with the boundary conditions $`W(0,t|x_0,0)=W(L,t|x_0,0)=0`$, and the initial condition $`W(x,0|x_0,0)=\delta (xx_0)`$ with $`0<\gamma 1`$. Equation (31) is straightforwardly solved by separation of variables :

$$W(x,t|x_0,0)=\frac{2}{L}\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{sin}\left(\frac{n\pi x_0}{L}\right)\mathrm{sin}\left(\frac{n\pi x}{L}\right)E_\gamma \left(K_\gamma \frac{n^2\pi ^2}{L^2}t^\gamma \right),$$
(32)
where we have taken into account that the solution of the fractional differential equation $`dT/dt=K_\gamma \lambda ^2{}_{0}{}^{}D_{t}^{1\gamma }T`$ is given by a Mittag-Leffler function $`T(t)=E_\gamma \left(K_\gamma \lambda ^2t^\gamma \right)`$. The survival probability of a subdiffusive random walker starting from $`x_0(0,L)`$ is given by $`_0^LW(x,t|x_0,0)𝑑x`$. Let $`\mathrm{\Phi }_L(t)`$ be the survival probability of the random walker averaged over all configurations of traps which contain the origin inside a hole of length $`L`$. Then $`\mathrm{\Phi }_L(t)=_0^L_0^LW(x,t|x_0,0)𝑑x𝑑x_0`$, and Eq. (32) yields

$$\mathrm{\Phi }_L(t)=\frac{8}{\pi ^2}\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{(2n+1)^2}E_\gamma \left(\frac{(2n+1)^2\pi ^2}{L^2}\lambda ^2K_\gamma t^\gamma \right).$$
(33)
The survival probability of the subdiffusive random walker in the two-sided trapping problem is finally given by an average of $`\mathrm{\Phi }_L(t)`$ over the distribution of hole lengths in the random trap configurations. This distribution is $`\eta (L)=\lambda ^2Le^{\lambda L}`$ , and consequently we have

$`\mathrm{\Phi }(t)`$ $`=\lambda ^2{\displaystyle _0^{\mathrm{}}}Le^{\lambda L}\mathrm{\Phi }_L(t)𝑑L`$
$`={\displaystyle \frac{8\lambda ^2}{\pi ^2}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(2n+1)^2}}{\displaystyle _0^{\mathrm{}}}Le^{\lambda L}E_\gamma \left({\displaystyle \frac{\beta _n}{L^2}}\right)𝑑L`$ (34)
with $`\beta _n=K_\gamma (2n+1)^2\pi ^2t^\gamma `$, $`n=1`$, $`2`$, $`\mathrm{}`$. This is the exact integral representation of the survival probability of the subdiffusive particle in the two-sided trapping problem.

### 4.1 Long-time behavior of $`\mathrm{\Phi }(t)`$

For $`\gamma =1`$ the Mittag-Leffler function reduces to an exponential and the integrals in Eq. (34) are estimated asymptotically for large times by means of the Laplace method . However, this approach is not possible for $`\gamma 1`$. In order to perform a long-time asymptotic evaluation of the integrals in Eq. (34), we use the following series expansion of the Mittag-Leffler function:

$$E_\gamma \left(\frac{\beta _n}{L^2}\right)=\underset{m=1}{\overset{p}{}}\frac{(1)^{m+1}L^{2m}}{\beta _n^m\mathrm{\Gamma }(1\gamma m)}+O\left(\frac{L^2}{\beta _n}\right)^{1+p},$$
(35)
which is valid for $`\beta _n/L^2\mathrm{}`$ in the range $`0<\gamma <2`$ . Inserting this series expansion into Eq. (34), integrating term by term, and after some rearrangements we find

$$\mathrm{\Phi }(t)=\frac{2}{\pi ^2}\underset{m=1}{\overset{p}{}}\frac{\left(42^{2m}\right)\zeta (2m+2)\mathrm{\Gamma }(2m+2)}{(1)^{m+1}\mathrm{\Gamma }(1\gamma m)\pi ^{2m}\xi ^{2m}}+O\left(\xi ^{2p2}\right),$$
(36)
where we have taken into account that $`_{n=0}^{\mathrm{}}\mathrm{\hspace{0.17em}1}/(2n+1)^{2m+2}=(12^{2m2})\zeta (2m+2)`$, $`\zeta (s)`$ being Riemann’s zeta function. The dominant term is

$$\mathrm{\Phi }(t)\frac{1}{2\mathrm{\Gamma }(1\gamma )\lambda ^2K_\gamma }t^\gamma ,$$
(37)
which is valid for $`\xi 1`$. Therefore we find an algebraic fluctuation slowdown of the survival probability corresponding to the Donsker-Varadhan  limit with an exponent two times greater than that of the one-sided long-time behavior in Eq. (14). The two-sided result $`\mathrm{\Phi }(t)t^\gamma `$ is well known and has been interpreted as an “avoided crossing effect”: trapping cannot be more efficient that the probability $`_t^{\mathrm{}}\psi (t)𝑑tt^\gamma `$ of remaining at the initial site (see  and references therein). However, note that this “avoided crossing effect” is absent for the one-sided trapping problem since the probability of trapping $`\mathrm{\Phi }(t)t^{\gamma /2}`$ decays even more slowly than the probability of remaining at the initial site.

### 4.2 Short-time behavior of $`\mathrm{\Phi }(t)`$

A simple analytical expression for the short-time behavior of $`\mathrm{\Phi }(t)`$ can also be derived. To do so, we calculate the time derivative of $`\mathrm{\Phi }_L(t)`$ in Eq. (33) obtaining

$$\frac{d}{dt}\mathrm{\Phi }_L(t)=\frac{8}{L^2}K_\gamma {}_{0}{}^{}D_{t}^{1\gamma }J(\gamma ),$$
(38)
where we have taken into account that Mittag-Leffler functions are the solution of the fractional relaxation equation  and

$$J(\gamma )=\underset{n=0}{\overset{\mathrm{}}{}}E_\gamma \left(4a^2(n+1/2)^2\right)$$
(39)
with $`a^2=K_\gamma \pi ^2t^\gamma /L^2`$. The major contribution to the sum defining $`J(\gamma )`$ in Eq. (39) when $`a`$ is small comes from large $`n`$. Hence, to find the lowest-order term in the small-$`a`$ behavior of $`J(\gamma )`$, we can approximate the sum by an integral:

$$J(\gamma )_0^{\mathrm{}}E_\gamma [4a^2x^2]𝑑x=\frac{1}{2a}\stackrel{~}{E}_\gamma (u=0),$$
(40)
where $`\stackrel{~}{E}_\gamma (u)=_0^{\mathrm{}}e^{uy}E_\gamma (y^2)𝑑y`$ is the Laplace transform of $`E_\gamma (y^2)`$. To calculate $`\stackrel{~}{E}_\gamma (u)`$ we will exploit the relation between the Mittag-Leffler and the Fox $`H`$ functions :

$$E_\gamma (y^2)=\frac{1}{2}H_{12}^{11}\left[y|\begin{array}{c}(0,1/2)\hfill \\ (0,1/2),(0,\gamma /2)\hfill \end{array}\right].$$
(41)
From Eq. (41) and the properties of Fox functions  we find

$$\stackrel{~}{E}_\gamma (s)=\frac{1}{2}H_{22}^{21}\left[s|\begin{array}{c}(1/2,1/2),(1\gamma /2,\gamma /2)\hfill \\ (0,1),(1/2,1/2)\hfill \end{array}\right].$$
(42)
Inserting the result in Eq. (42) into Eq. (40), we finally arrive at the following explicit expression for $`J(\gamma )`$:

$$J(\gamma )\frac{1}{4a}\frac{\pi }{\mathrm{\Gamma }(1\gamma /2)}=\frac{1}{4\mathrm{\Gamma }(1\gamma /2)}\frac{L}{\sqrt{K_\gamma }}t^{\gamma /2},$$
(43)
where we have calculated the value of the Laplace transform for $`s=0`$ by resorting to the general series expansion of Fox functions . From Eqs. (38) and (43) we can also write

$$\frac{d}{dt}\mathrm{\Phi }_L(t)\frac{2}{\mathrm{\Gamma }(\gamma /2)}\frac{\sqrt{K_\gamma }}{L}t^{\gamma /21},$$
(44)
where the fractional derivative has been performed according to the rule: $`{}_{0}{}^{}D_{t}^{\nu }t^\mu =\mathrm{\Gamma }(1+\mu )t^{\mu \nu }/\mathrm{\Gamma }(1+\mu \nu )`$ which is valid for arbitrary real parameters $`\mu `$ and $`\nu `$ . By integrating Eq. (44) with the initial condition $`\mathrm{\Phi }_L(t=0)=1`$ we find

$$\mathrm{\Phi }_L(t)1\frac{4}{\gamma \mathrm{\Gamma }(\gamma /2)}\frac{\sqrt{K_\gamma t^\gamma }}{L},$$
(45)
and performing the average in Eq. (34) over trap configurations we finally obtain the short-time behavior of the survival probability of a subdiffusive random walker in the two-sided trapping problem

$$\mathrm{\Phi }(t)1\frac{4\xi }{\gamma \mathrm{\Gamma }(\gamma /2)}.$$
(46)
or, for $`\xi 1`$,

$$\mathrm{\Phi }(t)\mathrm{exp}\left(\frac{2\xi }{\mathrm{\Gamma }(1+\gamma /2)}\right).$$
(47)
This expression is just the zeroth-order Rosenstock approximation $`\mathrm{\Phi }(t)=\mathrm{exp}[\lambda S_{}(t)]`$ \[see Eq. (16)\] because the average two-sided span $`S_{}(t)`$ is simply twice the average one-sided span $`S(t)`$ and, by Eq. (26), $`S(t)=(K_\gamma t^\gamma )^{1/2}/\mathrm{\Gamma }(1+\gamma /2)`$. In Fig. (2) we compare simulation results for subdiffusive random walkers in the case $`\gamma =1/2`$ with the numerical integration of the exact integral representation of $`\mathrm{\Phi }(t)`$ given by Eq. (34), the long-time behavior given in Eq. (36), and the Rosentock short-time behavior in Eq. (47).

## 5 CONCLUSIONS AND REMARKS

In this paper we considered the one-dimensional trapping problem for subdiffusive particles. In this problem a single subdiffusive random walker starting at $`x=0`$ on a lattice randomly filled with absorbing static traps with a density $`c`$ performs a random exploration until he encounters one of these traps. The quantity of main interest is the survival probability $`\varphi (t)`$ of the particle. We have considered two versions of the problem: the one-sided case in which only one half-line is filled with traps, and the two-sided case corresponding to a random filling of both sides of the line $`x>0`$ and $`x<0`$. In the context of ordinary diffusive random walkers, this problem has a long tradition, and many applications to physics and chemistry have been discussed .

Great interest has also arisen recently around subdiffusive anomalous processes (for example, as a way of mimicking transport in disordered media), and a continuous fractional diffusion description has been put forward (see  and references therein).It thus seems convenient to extend the fruitful trapping model to the case of subdiffusive particles. We achieved this objective for the trapping of a single subdiffusive random walker in one dimension and, by means of the fractional diffusion formalism, exact analytical expressions were found in terms of the special functions characteristic of fractional calculus. There is a possibility for these results to be checked in diffusion experiments performed in constrained geometries and/or disordered materials . From another point of view, the trapping problem could also be interpreted as a pseudo-first-order reaction of the form $`A+TT`$, where $`A`$ is the random walker and $`T`$ is the trap. There is also an increasing interest in reaction-subdiffusion processes , and our exact results apply to a special class of these processes.

The long-time behavior of the survival probability $`\mathrm{\Phi }(t)`$ has been a subject of particular interest in the ordinary diffusive trapping problem . In this regime an anomalous fluctuation slowdown, known as the Donsker-Varadhan limit, has been reported. In the one-sided and two-sided subdiffusive models we found the algebraic decays $`\mathrm{\Phi }(t)t^{\gamma /2}`$ and $`\mathrm{\Phi }(t)t^\gamma `$, respectively. For normal diffusive particles, the algebraic decay $`\mathrm{\Phi }(t)t^{1/2}`$ was found in  for the one-sided trapping case and interpreted as due to the long explorations carried by the random walkers on the half-line free of traps . A stretched exponential decay is the corresponding behaviour for the long-time two-sided trapping of diffusive particles . The two-sided subdiffusive result $`\mathrm{\Phi }(t)t^\gamma `$ is well known and has been interpreted as an “avoided crossing effect” (see  and references therein). This effect is absent for the one-sided trapping problem since the probability of trapping $`\mathrm{\Phi }(t)t^{\gamma /2}`$ for this case decays even more slowly than the probability of remaining at the initial site.

The present work can continue along several directions. First, the generalization of the subdiffusive trapping problem to the case of $`N>1`$ independent subdiffusive random walkers could be analyzed as has already been done in the diffusive case . Application of the same techniques to higher dimensional spaces is also an obvious approach, but we consider it unlikely to produce exact closed expressions for the relevant quantities. Trapping models of subdiffusive particles in which the trapping process is stochastically “gated” is also an interesting field. The effect of “gates” in diffusion-limited reactions has recently been considered in connection with reactions of complex molecules in biological media . The chemically active groups of these molecules may be screened by the inactive parts, giving rise to effective reactivities described by Poisson processes. The media where these reactions takes place are usually disordered and reaction-subdiffusion models should provide a better description of these processes.

This work was supported by the Ministerio de Ciencia y Tecnología (Spain) through Grant No. BFM2001-0718 and by the European Community’s Human Potential Programme under contract HPRN-CT-2002-00307, DYGLAGEMEM.