text:
# Measurement of the Sunyaev-Zel’dovich Increment in massive galaxy clusters

## 1 Introduction

The Sunyaev-Zel’dovich (SZ) effect can be used as a powerful probe of cosmology. Because the intensity of the SZ distortion is virtually independent of redshift, SZ measurements provide an avenue to study clusters and their peculiar velocities to any distance. Additionally, the SZ effect allows independent determination of various cosmological parameters. For example, because the SZ effect’s amplitude is proportional to the electron density along the line of sight through the cluster, $`n_\mathrm{e}`$, while the X-ray amplitude is proportional to $`n_\mathrm{e}^2`$, SZ effect data can yield distances to clusters independently from the cosmological distance ladder. This information can then be used to provide estimates of the Hubble constant $`H_0`$ (e.g. Reese et al. 2000).

The physics of the SZ effect is quite simple (Sunyaev & Zel’dovich 1972, Birkinshaw 1999, Carlstrom, Holder & Reese 2002). Cosmic microwave background (CMB) photons can scatter off hot ($`T_\mathrm{e}10^7`$K) electrons in a plasma. A statistical net gain of energy of the photons from the electrons occurs, producing a characteristic distortion in the CMB as seen through the plasma. This distortion produces a decrement in the CMB’s temperature below about 200 GHz, but an increment in the CMB’s temperature above this frequency.

The most important use of SZ increment measurements lies in the determination of the spectral shape of the SZ effect in clusters of galaxies. This shape is the sum of the thermal and kinetic effects together with other sources of emission within the cluster. The kinetic SZ effect is a change in the apparent CMB temperature as seen through the cluster due to the peculiar motion of the cluster relative to the CMB. The kinetic SZ effect may be present in any given cluster and may manifest itself as either a temperature increase or decrease, although it never dominates the thermal effect for expected cluster velocities (Birkinshaw, 1999). Knowledge of the cluster’s properties obtained at several wavelengths can allow separation of these two effects. If the kinetic effect can be isolated, the cluster’s peculiar velocity along the line of sight can be determined. Statistics about large scale flows at a variety of redshifts provide strong constraints on the dynamics of structure formation. In principle, detailed measurement of the SZ spectral shape can also constrain the cluster gas temperature through the relativistic correction (Rephaeli 1995, Birkinshaw 1999, Colafrancesco, Marchegiani, & Palladino 2003).

High frequency measurements can also help separate contaminants such as dusty galaxies which may exist in millimetric measurements, since the dominant sub-mm point sources differ from those seen at longer wavelengths.

Measurement of the Sunyaev-Zel’dovich (SZ) effect increment in galaxy clusters has been an elusive and important goal of observational cosmology for over a decade. While measurements of the SZ effect decrement are common (Birkinshaw 1999, Carlstrom, Holder & Reese 2002), measurement and study of the increment has been less successful. This is largely because current instruments generally lack the sensitivity to measure small amplitude, extended, positive signals in the sky. Also, other positive flux sources in a cluster field can confuse SZ measurements (Loeb & Refregier 1997, Blain 1998). However, experiments working near the null of the thermal SZ effect like SuZIE (Holzapfel et al., 1997), and ACBAR (Peterson et al., 2002) are beginning to yield results. The balloon-borne experiment PRONAOS has had limited success in SZ increment measurements (Lamarre et al., 1998), and other detections are also claimed (e.g. Komatsu et al. 1999). Recently, Benson et al. (2003) have performed this type of measurement in six galaxy clusters, and find that measurement of any given cluster’s peculiar velocity is difficult.

Measurement of the SZ increment is difficult for a variety of reasons, including the amplitude of the atmospheric signal at sub-mm wavelengths. $`450`$$`\mu `$m and $`850`$$`\mu `$m observations take advantage of windows in the opacity of the sky at these wavelengths, but successful observation of any astronomical source in the sub-mm regime still requires non-trivial atmospheric subtraction techniques. Additionally, because the SZ effect is a small amplitude signal extended across a relatively large region of the sky, sensitive instruments which can sample a range of spatial frequencies are required. Care must be taken when reconstructing the spatial shape of the SZ effect from data recorded as differences between two or more telescope pointings. As point sources are expected to be a major contaminant in the sub-mm regime, high resolution is also required to adequately understand their effect in a given field. One of the instruments currently most capable of performing observations of the SZ increment is the Sub-millimetre Common User Bolometer Array (SCUBA; Holland et al. 1999). This instrument provides enough sensitivity and angular resolution to detect and remove point sources (including lensed background sources) from the data. SCUBA also samples large enough regions of the sky to constrain the shape of the SZ increment in moderate to high redshift clusters. However, even with these characteristics, measurement of the SZ increment with SCUBA is not trivial.

## 2 Observations

The bulk of the data discussed here was obtained on November $`10`$ and $`11`$, $`1999`$ at the James Clerk Maxwell Telescope (JCMT) at the summit of Mauna Kea, Hawaii. The JCMT combines a high, dry site with a $`15`$ metre dish and sensitive instruments. SCUBA, attached on the Nasmyth platform of the JCMT, allows simultaneous observations with $`91`$ $`450`$$`\mu `$m–band bolometers and $`37`$ $`850`$$`\mu `$m–band bolometers. Because the $`450`$$`\mu `$m band observes very little SZ increment flux, while the expected peak of the SZ signal occurs near $`850`$$`\mu `$m, the $`450`$$`\mu `$m data can act as a monitor of atmospheric emission in this experiment.

SCUBA is sparsely sampled spatially, which is undesirable in this experiment. ‘Jiggling’ to fully sample the array consists of small ($`6`$ arcsec) changes in the position of the beam in a regular pattern. The detectors integrate in each jiggle position for $`1.024`$s. Although $`64`$ positions are required to fully sample both arrays, this experiment used $`16`$ positions because it was only necessary to fully sample the $`850`$$`\mu `$m array. Moreover, executing the jiggle pattern as quickly as possible is more conservative in terms of residual sky fluctuations.

Variations in the signal produced by atmospheric (sky) noise and thermal offsets from the instrument are removed by rapid chopping together with slower telescope nodding (e.g. Archibald et al. 2002). Rapidly chopping the secondary mirror faster than the rate at which the sky is varying removes much of this atmospheric signal. Nodding the telescope’s primary mirror to a second reference position removes slow sky noise gradients and instrument noise. These motions produce a chop pattern such that the two reference beams are each observed half as often as the source beam (Fig. 1). In the standard SCUBA mapping mode, the chop occurs at about $`8`$ Hz and the nod occurs at about $`1/16`$ Hz. Although a chop throw in the range $`50`$ to $`100`$ arcsec is standard, this experiment requires a chop that does not put significant SZ flux in the reference beams (Fig. 1). For this reason, an extremely large JCMT secondary mirror chop size of $`180`$ arcsec was used.

$`89`$ks of data were taken, split between several fields. The CSO $`225`$GHz $`\tau `$-meter was not operational for either observing shift, so we relied on skydips to measure the atmospheric optical depth, $`\tau _{850}`$, in our data set. Both nights had relatively constant $`\tau _{850}`$, but the first evening had a lower average ($`\tau _{850}0.15`$) than the second ($`\tau _{850}0.3`$). Uranus was used as a calibration source throughout.

Three clusters and three blank fields were observed during this run (see Table 1 for details), referred to as Cl$`\mathrm{\hspace{0.17em}0016}`$, MS$`\mathrm{\hspace{0.17em}0451}`$, MS$`\mathrm{\hspace{0.17em}1054}`$, Blank $`02`$, Blank $`09`$ and Blank $`25`$ hereafter. $`23`$ks of additional data were collected with an aluminum reflector masking the aperture of SCUBA. These data are used as null tests to study SCUBA’s behaviour under both optical and electrical influences (blank sky data) and electrical influences only (reflector data).

Additional data were taken at the JCMT on October $`2`$, $`3`$, and $`6`$, $`2002`$ in the same manner as the $`1999`$ run. Both Cl$`\mathrm{\hspace{0.17em}0016}`$ and another blank field, Blank $`21`$, were mapped. In addition to the jiggle mapping, photometry of possible point sources in the Cl$`\mathrm{\hspace{0.17em}0016}`$ field were performed. The photometry used a $`9`$ point jiggle pattern with $`1`$ arcsec spacing and a $`60`$ arcsec chop throw. Bolometer G9 of the SCUBA array was not operating properly during this run and its data are not used here. A total of $`13`$ks of data were taken in 2002; how they are combined with the earlier data is discussed below.

## 3 General Analysis

### 3.1 Data Preparation

The data are deswitched, flat-fielded and corrected for atmospheric extinction with surf (the SCUBA User Reduction Facility, Jenness & Lightfoot 1998). Extinction correction is performed by calculating the airmass at which each bolometer measurement was made, and then multiplying by the zenith sky extinction at the time of the measurement to give the extinction optical depth along the line of sight. Each data point is then multiplied by the exponential of the optical depth to give the value that would have been measured in the absence of the atmosphere. The photometric calibration of deep SCUBA fields like those presented here is known to be accurate to about 15 per cent (Borys 2002, Archibald et al. 2002), although $`2/3`$ of this uncertainty is due to error in the absolute calibration of the planets’ flux.

It is well documented that the deswitched $`450`$$`\mu `$m and $`850`$$`\mu `$m time streams are highly correlated, and that this correlation stems primarily from atmospheric noise (e.g. Jenness, Lightfoot, & Holland 1998, Borys et al. 1999, Archibald et al. 2002). The standard method of removing atmospheric noise from an array is to subtract the array average from each bolometer at each time step. However, most of the information about the SZ flux in this measurement is contained in the difference between the source and reference beam fluxes (see Fig. $`1`$) rather than in the shape of the SZ emission across the field of view. Therefore, removing the $`850`$$`\mu `$m array average at each time step would drastically reduce the inferred SZ emission. However, because the SZ effect is small at $`450`$$`\mu `$m, using that band’s data should mitigate this problem. A least-squares linear fit of the $`450`$$`\mu `$m array average to the $`850`$$`\mu `$m array average at each time step is used to obtain a linear coefficient $`C_1`$ and an offset $`C_0`$ (hereafter referred to as the $`450\mu \mathrm{m}`$ method). These fit parameters and the $`450`$$`\mu `$m array averages are then used in place of the $`850`$$`\mu `$m array average to subtract the atmospheric noise at each time:

$$S_b^{850}(t)=R_b^{850}(t)\left(C_1R^{450}(t)+C_0\right),$$
(1)
where $`S_b^{850}(t)`$ is the corrected data at time step $`t`$ for bolometer $`b`$, $`R_b^{850}(t)`$ is the raw, double-differenced and extinction corrected data at time step $`t`$ for bolometer $`b`$, the angled brackets denote an average over the whole array, and the superscripts denote wavelength. The fits are performed separately for each (approximately $`20`$ minute) observation. $`C_1`$ removes the atmospheric emission on scales of about $`1`$s, while $`C_0`$ removes the atmospheric average over periods long compared to the telescope nod time. The value of $`C_1`$ and $`C_0`$ do not vary by more than 10 per cent over an observation; this is consistent with the results of Borys et al. (1999). $`C_1`$ is typically $`0.6`$ and $`C_0`$ is of order $`1\mu `$V in our data set. The linear Pearson correlation coefficients of the $`450`$$`\mu `$m and the $`850`$$`\mu `$m time streams are on average $`0.95`$ and never less than $`0.90`$ for these observations. The average ratio of the variance of $`S_b^{850}`$ in the time stream of the $`450\mu \mathrm{m}`$ atmospheric subtraction method to the variance of the $`850`$$`\mu `$m time stream in the standard surf atmospheric subtraction method is $`1.15`$. This implies that the $`450\mu \mathrm{m}`$ method is slightly inferior for removing atmospheric noise. However, this approach is required for the analysis of SZ data, since it allows retention of information on angular scales up to that of the chop throw.

Of the seven target fields, three are unusable for this experiment in some way. A fundamental requirement of this analysis method is high enough signal to noise ratio that it is possible to determine the location of point sources in a field and remove them. However, the Blank $`02`$ field had only a short integration time (see Table 1), hence it is not included in the results of this experiment. Blank $`22`$ was a poor choice for a control field, because it contains large point sources which corrupt it for our purposes; Chapman et al. (2001) examine this field in detail. The MS$`\mathrm{\hspace{0.17em}0451}`$ field contains a sub-mm bright gravitationally lensed arc serendipitously discovered when combining this data with previous observations (Chapman et al., 2002). Although an interesting target in its own right, the $`12`$mJy beam<sup>-1</sup> resolved arc dominates the sub-mm emission in this cluster (Borys et al., 2003), making measurement of the SZ increment impossible with these data.

### 3.2 Point Source Removal

Recognition and removal of point sources, especially unbiased handling of lensed background sources, is important in this experiment (Loeb & Refregier 1997, Blain 1998). We adopt the following procedure. The measured flux differences are binned according to the spatial position of the ‘source’ beam. The standard method of presenting a SCUBA image is to convolve the binned data set with the telescope point spread function (PSF). Maps made by this method are shown in Fig. $`2`$ in order to indicate the locations of possible point sources and to exhibit their consistency with the same fields in Chapman et al. (2002). However, these images are not used as part of the quantitative analysis.

Instead, the binned flux differences are convolved with a wavelet similar to a ‘Mexican Hat’ function:

$$\psi (\theta )=\psi _0\left(1\frac{\theta ^2}{2\sigma _\mathrm{w}^2}\right)\mathrm{exp}\left(\frac{\theta ^2}{2\sigma _\mathrm{w}^2}\right).$$
(2)
Here $`\theta `$ is an angular distance and $`\sigma _\mathrm{w}`$ is a parameter charactering the wavelet’s width. The normalization $`\psi _0`$ is set such that the flux from a point source placed in a blank field is returned identically at $`\theta =0`$. This wavelet is more useful than the telescope PSF in this analysis because it returns an amplitude which is insensitive to the average value in the map, which is not expected to be zero in an SZ experiment. The wavelet’s width was chosen to be $`\sigma _\mathrm{w}=6.3`$ arcsec, slightly smaller than the $`\sigma `$ corresponding to the FWHM of the telescope, because this value maximizes the signal to noise ratio in simulated measurements. After convolution with this function, the positive and negative extrema in the maps are the locations of possible point sources in the source and reference beams respectively (Table 2).

The philosophy of this experiment’s point source removal scheme is not to develop a list of point sources per se, but rather to remove the flux contributed by any point-like source. In the simple method utilized here, possible point sources in a field are removed by adopting a single flux level above which sources are found. The limiting flux level is estimated based on atmospheric extinction, integration time and the noise equivalent flux density. For each field, we use the following thresholds: 9 mJy in Cl$`\mathrm{\hspace{0.33em}0016}`$, 9 mJy in MS$`\mathrm{\hspace{0.17em}1054}`$; 11.5 mJy in Blank $`09`$; and 13 mJy in Blank $`21`$. It is found that the point source list matches that found in Chapman et al. (2002) well. The estimated variance given by this method and the actual mean variance of the pixels agree to 10 per cent or better in each field.

By determining the reference beam positions at which most of the negative flux originates, it is found that the ‘two’ negative point sources in MS$`\mathrm{\hspace{0.17em}1054}`$ are, in fact, probably caused by one bright source. The source appears to be two sources because the MS$`\mathrm{\hspace{0.17em}1054}`$ field was observed only at the beginning and the end of the Nov. 11 1999 shift. The field rotated during the interval, thereby extending the single source into ‘two’ apparent sources. Our best estimate of this source’s true position is 10:57:03.5 $`\pm 2.0`$ seconds in right ascension and $`3`$:39:00 $`\pm 60`$ arc seconds in declination. Chapman et al. (2002) find another source of similar amplitude in this cluster which did not happen to fall in any of our beams. The presence of these sources is consistent with the source count models (e.g. Smail et al. 2002).

The effect of greater than 3$`\sigma `$ positive sources is removed directly from the time series using the measured amplitude, the telescope pointing history and the telescope’s PSF. A negative 3$`\sigma `$ source is equally likely to be in either of the two reference beam positions, so the same approach is not available. Hence, data associated with detected negative sources are simply excised from further analysis.

Monte Carlo simulations have been performed to develop an understanding of the effects of point sources on this data. These simulations are fully explained in Zemcov, Newbury & Halpern (2003) (hereafter ZNH ), and the application of the results found there are discussed in Sections $`3.3`$ and $`5`$.

### 3.3 Model Fitting

Because the SZ effect is extended on the sky, the simplest way to isolate it in a cluster is to fit the data to a model. Since the signal to noise ratio in any given pixel is poor for this data set, a simple isothermal $`\beta `$ model is adequate:

$$\mathrm{\Delta }I^\mathrm{D}(\theta )=\mathrm{\Delta }I_0^\mathrm{D}\left(1+\frac{\theta ^2}{\theta _0^2}\right)^{(13\beta )/2}.$$
(3)
Here, $`\mathrm{\Delta }I_0^\mathrm{D}`$ is a measure of the magnitude of the increment at the centre of the cluster, $`\theta `$ is the angular distance from the centre, and $`\theta _0`$ and $`\beta `$ are parameters characterizing the cluster. These parameters have already been determined for these clusters from combinations of X-ray and SZ decrement data (see Table $`3`$). Although there is some uncertainty in the values of $`\beta `$ and $`\theta _0`$, the error this causes in the inferred $`\mathrm{\Delta }I_0^\mathrm{D}`$ is negligible at this experiment’s level of precision. We thus fix the angular dependence and focus on determining the amplitude alone with our SCUBA data.

Because the SZ increment is non-zero at $`450`$$`\mu `$m and the $`450`$$`\mu `$m array average is used to remove residual sky noise, the superscript ‘D’ is used as a reminder that the value quoted is actually a difference between the $`850`$$`\mu `$m signal and the $`450`$$`\mu `$m array average. Retrieving the $`850`$$`\mu `$m signal (denoted $`\mathrm{\Delta }I_0^{850}`$) from $`\mathrm{\Delta }I_0^\mathrm{D}`$ will be discussed further below.

Directly fitting the data to an isothermal $`\beta `$ model differenced according to the telescope pointing history is the most straightforward way to determine the peak SZ increment flux in a given field. Both the $`\mathrm{\Delta }I_0^\mathrm{D}`$ and the 1$`\sigma `$ error in $`\mathrm{\Delta }I_0^\mathrm{D}`$ found from fitting to the differenced isothermal $`\beta `$ model are listed in Table $`4`$. The error in $`\mathrm{\Delta }I_0^\mathrm{D}`$ is found via the correlation matrix of the linear fit. $`\sigma (\mathrm{\Delta }I_0^\mathrm{D})`$ is the square root of the second diagonal element of this correlation matrix. The fit correlation matrix has very small off diagonal elements, meaning that the fit is not correlated with the (meaningless) voltage offset value.

To test the hypothesis that much of the information about the SZ intensity is contained in the difference $`S_b^{850}(t)`$, its average is calculated over the array at each time step. The mean of these values should be greater than 0 and be detected with approximately the same or slightly less significance as the detection given by the fit to the isothermal $`\beta `$ model. It is found that the means of the Cl$`\mathrm{\hspace{0.17em}0016}`$ and MS$`\mathrm{\hspace{0.17em}1054}`$ fields are approximately $`1.5\sigma `$ and $`2\sigma `$ above zero, respectively. We conclude that merely using the shape of the SZ distortion across the SCUBA field of view underestimates its true amplitude.

Note that intensities, $`\mathrm{\Delta }I_{\mathrm{SZE}}`$, are generally quoted in mJy per JCMT beamsize, which is related to temperature change via the derivative of the blackbody distribution:

$$\mathrm{\Delta }I_{\mathrm{SZE}}=2k_\mathrm{B}^3\mathrm{\Omega }\left(\frac{T_{\mathrm{CMB}}}{hc}\right)^2\frac{x^4e^x}{(e^x1)^2}\mathrm{\Delta }T_{\mathrm{SZE}}$$
(4)
where $`x=h\nu /kT_{\mathrm{CMB}}`$ is the dimensionless frequency, $`T_{\mathrm{CMB}}=2.725`$K, $`\mathrm{\Delta }T_{\mathrm{SZE}}`$ is the change in CMB temperature due to the SZ effect, and $`\mathrm{\Omega }=5.74\times 10^9`$ steradians per beamsize.

It is also desirable to make a crude ‘map’ to check the radial profile of the data. A matrix inversion technique utilizes all of the available data to make a map, including information about the off-source pointings (e.g. Tegmark & Bunn 1995, Wright, Hinshaw & Bennett 1996, Stompor et al 2002, Hinshaw et al. 2003). This method involves dividing the sky into a vector of pixels, creating matrices based on pointing and flux data, and the inversion of a single matrix to create a map. The matrix inversion itself is performed using a Singular Value Decomposition (e.g. Press et al. 1992).

Here, the pixel model of the sky is chosen to be the brightnesses of annular rings centred on the cluster, as this conforms to the approximately circular symmetry expected in the SZ signal. The rings have $`40`$ arcsec widths to track the radial variation with reasonable signal to noise ratio. Once the fluxes of the annuli including their errors and correlations have been estimated, they can be fit to the isothermal $`\beta `$ model. The fit is performed using both data and model with a weighted mean of 0. This ensures that the model’s DC term does not falsely increase the measured increment. A generalized $`\chi ^2`$ statistic is determined, where $`𝐱_\mathrm{m}`$ is a pixel’s isothermal $`\beta `$ model value, and $`𝐱_\mathrm{d}`$ is a pixel’s data value.

$$\chi ^2=(𝐱_\mathrm{m}𝐱_\mathrm{d})_i𝐂_{ij}^1(𝐱_\mathrm{m}𝐱_\mathrm{d})_j.$$
(5)
Our observational strategy produces non-negligible correlations between the annular pixels which are given in $`𝐂_{ij}`$, the pixel–pixel covariance matrix. Explicitly, $`𝐂_{ij}=\left(𝐏^\mathrm{T}𝐍^1𝐏\right)^1`$, where $`𝐏`$ is our pointing matrix and the noise is $`𝐍`$. $`𝐍`$ is assumed to be diagonal because removing the fit to the $`450`$$`\mu `$m array average effectively acts as a pre-whitening filter. These correlations need to be carefully included when the SZ amplitude is determined using this annular pixel model and an isothermal $`\beta `$ profile. Fig. $`3`$ shows the annular pixel fluxes as well as the best fit profile (Equation ($`3`$)) using the parameters listed in Tables $`3`$ and $`4`$. The two methods yield very similar results, although we find that the direct fitting method is more sensitive.

Generally, the figure of merit for the goodness of fit is that $`\chi ^2`$ should be approximately the number of degrees of freedom in the data. However, source confusion may have an effect on the applicability of this figure of merit. To test how well $`\chi ^2`$ estimates the goodness of fit, Monte Carlo simulations similar to those presented in ZNH  are performed in order to determine the probability that the cluster field fits are consistent with the absence of an SZ increment. The direct fit to the isothermal $`\beta `$ model differences method is used to find $`\mathrm{\Delta }I_0^\mathrm{D}`$ in simulations containing realistic point sources and with noise levels set to be the same as those in this experiment, but with no SZ increment in the fields. Point sources are removed at the same level as in this experiment; Fig. $`4`$ shows the results of 200 such simulations. It is found that the fits associated with the highest $`10`$ per cent of the $`\chi ^2`$ statistics always yield drastically incorrect increment values because of point source contamination. We therefore exclude the worst $`10`$ per cent of the fits. Fig. $`4`$ also shows the SCUBA results for both clusters, each of which pass the $`\chi ^2`$ cutoff. The results overlap, as the fits for these fields give almost the same increment value. If the SZ increment were absent from Cl$`\mathrm{\hspace{0.17em}0016}`$ and MS$`\mathrm{\hspace{0.17em}1054}`$, the probability of obtaining the $`\mathrm{\Delta }I_0^\mathrm{D}`$ found is less than approximately $`0.01`$ for either cluster.

The $`\chi ^2`$ values found for Cl$`\mathrm{\hspace{0.17em}0016}`$ are $`3.7`$ for $`6`$ degrees of freedom via the seven pixel fit method, and a $`\chi ^2`$ of $`178402`$ for $`178406`$ degrees of freedom using the direct fit to model differences. We conclude that the fits are fully consistent with the presence of an increment in this cluster.

Our confidence in the parameters found for MS$`\mathrm{\hspace{0.17em}1054}`$ is less than that for Cl$`\mathrm{\hspace{0.17em}0016}`$. Based on its $`\chi ^2`$ of $`13.4`$ for $`6`$ degrees of freedom from fitting the annular pixels to an isothermal $`\beta `$ model, this cluster would normally be flagged as a poor fit. The fit may be poor because the isothermal $`\beta `$ model does not describe the distribution of MS$`\mathrm{\hspace{0.17em}1054}`$’s intracluster electron gas particularly well (Jeltema et al., 2001). However, the fit may become tolerable if the nature of the double–cored electron gas distribution is taken into account. Nevertheless, the reduced $`\chi ^2`$ from fitting to the isothermal $`\beta `$ model differences is $`258438`$ for $`258204`$ degrees of freedom in this cluster. The reduced $`\chi ^2`$ statistic associated with this fit falls well below the worst 10 per cent of reduced $`\chi ^2`$ statistics in the Monte Carlo simulations discussed above. The fit is therefore acceptable based on this criterion.

The parameters determined for the blank fields are also given in Table $`4`$. From the fits to annular pixels, the Blank $`09`$ field has a $`\chi ^2`$ of $`11.7`$ for $`6`$ degrees of freedom and the Blank $`21`$ field has a $`\chi ^2`$ of $`13.6`$ for $`6`$ degrees of freedom. The $`\chi ^2`$ from fitting to the model differences are poor in both fields. Both fields also seem to exhibit unremoved, confused point sources. The simulations of ZNH  show that fields with no SZ increment form a well-defined locus on the $`I_0`$$`\chi ^2`$ plane, and both of the blank fields in this experiment would be inside the confused regime and hence be rejected as SZ increment candidates based on their $`\chi ^2`$.

For the purposes of comparison, the line $`\mathrm{\Delta }I_0=0`$mJy beam<sup>-1</sup> is fit to the data using the direct fit method. The $`\mathrm{\Delta }\chi ^2(=\chi ^2\chi _{\mathrm{min}}^2)`$ values for these fits are 6.7 for Cl$`\mathrm{\hspace{0.17em}0016}`$, 6.6 for MS$`\mathrm{\hspace{0.17em}1054}`$, 0.5 for Blank $`09`$ and 2.0 for Blank $`21`$. These values are in good agreement with the significance of the detections.

Although the MS$`\mathrm{\hspace{0.17em}0451}`$ and Blank $`22`$ fields were found to be unsatisfactory for this experiment earlier in our analysis, this analysis pipeline also rejects these fields. Both fields have $`\mathrm{\Delta }I_0`$ consistent with $`0`$mJy per beam, and poor $`\chi ^2`$ due to the number and location of sources within them. This is an objective check that retrieving an SZ increment from data heavily polluted by point sources is difficult.

## 4 Individual Cluster Analysis & Results

In this section our data are used to determine the $`850`$$`\mu `$m SZ increment values in Cl$`\mathrm{\hspace{0.17em}0016}`$ and MS$`\mathrm{\hspace{0.17em}1054}`$. A maximum likelihood method is used to determine the best fit Compton $`y`$ parameter (which parametrizes the integrated pressure in the cluster along the line of sight) and, in Cl$`\mathrm{\hspace{0.17em}0016}`$, the peculiar velocity of the cluster along the line of sight.

### 4.1 Cl$`\mathbf{\hspace{0.17em}0016}`$+$`\mathrm{𝟏𝟔}`$ $`𝚫𝐈_\mathrm{𝟎}^{\mathrm{𝟖𝟓𝟎}}`$ Determination

Worrall & Birkinshaw (2003) find that Cl$`\mathrm{\hspace{0.17em}0016}`$ has $`T_\mathrm{e}9`$keV, implying that it has a relativistic electron population which may give a significant correction near the positive peak of the SZ spectrum. Therefore, the corrections found by Itoh, Kohyama & Nozawa (1998), which are good to fifth order in $`k_\mathrm{B}T_\mathrm{e}/m_\mathrm{e}c^2`$, are used to determine the correct shape of the thermal SZ distortion. These corrections, which represent a $`10`$$`20`$ per cent change to the flux at $`850`$$`\mu `$m for standard cluster parameters, are applied to all of the spectral calculations, and are implicit in the use of the term ‘thermal’ SZ effect below.

Because the SCUBA filters have a finite bandwidth, we need to determine the effective central frequencies of this measurement. This calculation must be performed because the SZ spectrum differs from the grey-body spectra SCUBA usually observes. To do this, the ‘equivalent frequency’ is defined as:

$$\overline{\nu }=\frac{_0^{\mathrm{}}\nu F(\nu )S(\nu )𝑑\nu }{_0^{\mathrm{}}F(\nu )S(\nu )𝑑\nu },$$
(6)
where $`F(\nu )`$ is the actual SCUBA filter’s response and $`S(\nu )`$ is the SZ effect spectrum over the region of interest. The results of this calculation are listed in Table $`5`$; we conclude that the frequency shift from the nominal values are negligibly small.

Foreground dust emission could also contaminate the measurement of the SZ increment in these clusters. ISO<sup>1</sup><sup>1</sup>1We have used archival data from the Infrared Space Observatory (ISO), a European Space Agency mission with the participation of ISAS and NASA (Kessler et al., 1996). data are used to estimate the dust contribution. The telescope pointings are used to sample the $`180`$$`\mu `$m and $`90`$$`\mu `$ISO map of the cluster region in exactly the same way as the SCUBA data. Errors for each measurement are determined from the error map and incorporated into the simulated time stream. These data are then fit to the model differences. The $`\chi ^2`$ values imply that the fits are very poor; this is expected because any dust emission in the region (Galactic cirrus, extragalactic point sources, or the cluster itself) is unlikely to be distributed like the isothermal $`\beta `$ model. However, the dust may still contribute to the SZ signal and so must be removed. The dust contribution is determined by fitting a modified black body of the form:

$$I_\mathrm{d}(\nu ,T_\mathrm{d})=\alpha \frac{\nu ^{3+\gamma }}{e^{h\nu /k_\mathrm{B}T_\mathrm{d}}1}.$$
(7)
Here, the parameters are the spectral index $`\gamma `$, dust temperature $`T_\mathrm{d}`$, and normalization $`\alpha `$, which is determined using both of the ISO data points weighted by their errors. $`\gamma `$ is taken to be $`2.0`$ and $`T_\mathrm{d}`$ to be $`20`$K, parameters which are typical for ISO–detected cirrus at these wavelengths (e.g. Lamarre et al. 1998). Using this model, the dust emission at $`450`$$`\mu `$m is determined to be $`I_\mathrm{d}(660\mathrm{GHz})=0.07\pm 0.04`$ mJy beam<sup>-1</sup> for Cl$`\mathrm{\hspace{0.17em}0016}`$. This value is not very sensitive to variations in the parameters $`\gamma `$ and $`T_\mathrm{d}`$. The fit is corrected for dust emission by the addition of this value, which is a few per cent correction.

Because our analysis scheme subtracts the $`450`$$`\mu `$m flux, a method of retrieving $`\mathrm{\Delta }I_0^{850}`$ from $`\mathrm{\Delta }I_0^\mathrm{D}`$ is required. The actual quantity subtracted from each $`850`$$`\mu `$m bolometer’s time stream is the average of the $`450`$$`\mu `$m array at each time step. This average contains the SZ signal sampled as described in Fig. 1, with $`\mathrm{\Delta }I_0^{450}`$ being the peak isothermal $`\beta `$ model amplitude. The average value of the $`450`$$`\mu `$m array is determined by sampling the model SZ profile in the same way as in this experiment. A map is made based on this data, and the average of this map in the field of view is constructed. It is found that there is a linear relationship between the input SZ amplitude and the average of the differenced map, given by $`\mathrm{\Delta }I_{\mathrm{Ave}}^{450}=A\mathrm{\Delta }I_0^{450}`$, where $`A`$ is a linear coefficient and $`\mathrm{\Delta }I_{\mathrm{Ave}}^{450}`$ is the average of the difference map over the field of view. This means that $`\mathrm{\Delta }I_0^{850}`$ can be found via:

$$\mathrm{\Delta }I_0^{850}=\mathrm{\Delta }I_0^D+A\mathrm{\Delta }I_0^{450},$$
(8)
where $`A`$ depends only on the shape of the isothermal $`\beta `$ model. $`A=0.72`$ in Cl$`\mathrm{\hspace{0.17em}0016}`$, including the gain difference between the $`850`$$`\mu `$m and $`450`$$`\mu `$m arrays.

To determine the best fit Compton $`y`$ parameter and peculiar velocity $`v_{\mathrm{pec}}`$, we utilize a likelihood function which searches the $`y`$ \- $`v_{\mathrm{pec}}`$ parameter space for its maximum. Previous decrement measurements at $`15`$GHz (Grainge et al., 2002), $`20`$GHz (Hughes & Birkinshaw, 1998) and $`30`$GHz (Reese et al., 2000) are used alone to determine a best-fit $`y`$ value, and also in conjunction with our measurement to help constrain both the thermal and kinetic SZ effects in Cl$`\mathrm{\hspace{0.17em}0016}`$. Alone, these values yield $`y=(2.38\pm 0.12)\times 10^4`$ in this cluster, assuming $`v_{\mathrm{pec}}=0`$km s<sup>-1</sup>. Our value alone yields $`y=(2.2\pm 0.7)\times 10^4`$ under the same assumption. This measurement is therefore consistent with earlier results.

We find that the $`850`$$`\mu `$m SZ increment value for Cl$`\mathrm{\hspace{0.17em}0016}`$ is $`\mathrm{\Delta }I_0^{850}=2.2\pm 0.7`$ mJy beam<sup>-1</sup>. This corresponds to $`0.38\pm 0.12`$ MJy sr<sup>-1</sup> or $`\mathrm{\Delta }T_{\mathrm{CMB}}=1.2\pm 0.4`$ mK in thermodynamic temperature units. Fig. $`5`$ shows the corrected $`850`$$`\mu `$m data point, along with the radio and ISO results. Also plotted are the best fit thermal SZ effect and dust emission spectra, and the position of the $`450`$$`\mu `$m band–centre for reference.

In addition to constraining $`y`$, the amplitude of the kinetic SZ effect is also constrained by this likelihood analysis. Using both our measurement and the three decrement measurements, it is found that the marginalized best fit parameters are $`y=(2.38\pm _{0.34}^{0.36})\times 10^4`$ and $`v_{\mathrm{pec}}=400\pm _{1400}^{1900}`$ km s<sup>-1</sup> in Cl$`\mathrm{\hspace{0.17em}0016}`$. Assuming $`v_{\mathrm{pec}}=0`$, these points yield $`y=(2.38\pm 0.11)\times 10^4`$. Fig. $`5`$ also shows the inferred kinetic SZ effect from the $`850`$$`\mu `$m data. This result is therefore only a weak constraint on the line of sight peculiar velocity.

A kinetic SZ effect in Cl$`\mathrm{\hspace{0.17em}0016}`$ may be mimicked by a primary fluctuation in the CMB, because the spectral distributions of the two effects are the same. However, the primary CMB anisotropies at the scales of this experiment, $`\mathrm{}>5000`$, are about $`1\mu `$K rms (Hu & Dodelson 2002, Borys et al. 1999). This corresponds to a $`\pm 0.005`$mJy beam<sup>-1</sup> signal, which is negligible at our level of precision.

### 4.2 MS$`\mathbf{\hspace{0.17em}1054}\mathrm{𝟎𝟑}`$ $`\mathrm{\Delta }I_0^{850}`$ Determination

The value of $`y`$ for MS$`\mathrm{\hspace{0.17em}1054}`$ is found by determining the difference $`I^{850}AI^{450}`$ for the SZ distortions defined by $`1.0\times 10^5y1.0\times 10^3`$, where $`A=0.81`$ for this cluster. The $`y`$ values giving differences most closely matching the data and 1$`\sigma `$ error points are the best fitting $`y`$ parameter and error bars. This method yields $`y=\left(2.0\pm 1.0\right)\times 10^4`$ for this cluster, including relativistic corrections. These $`y`$ values correspond to an increment of $`\mathrm{\Delta }I_0^{850}=2.0\pm 1.0`$mJy beam<sup>-1</sup>, or $`\mathrm{\Delta }T_{\mathrm{CMB}}=1.0\pm 0.5`$mK.

Because little other relevant data exist for this cluster, further analysis is not possible. The dust contamination cannot be estimated because there are no ISO data, and the available dust maps (Schlegel, Finkbeiner & Davis, 1998), which are of poor angular resolution, provide little spectral information. However, it appears that the dust contamination is at a similar level to that found in Cl$`\mathrm{\hspace{0.17em}0016}`$, and hence is probably negligible. Although one published estimate of the thermal decrement magnitude already exists for MS$`\mathrm{\hspace{0.17em}1054}`$ (Joy et al., 2001), a worthwhile estimate of its kinetic effect cannot be made because no explicit $`y`$ parameter or $`\mathrm{\Delta }T_{\mathrm{CMB}}`$ is available.

## 5 Discussion

In principle, measurements of the SZ increment, in combination with lower frequency measurements, can yield a great deal of information. For example, the line of sight peculiar velocity and possibly even the cluster temperature can be determined via measurement of the SZ increment. However, measurement of the SZ increment with SCUBA is difficult. There are three essential components which combined to make this experiment a success.

The first component is a carefully planned observational strategy designed to maximize SCUBA’s sensitivity to an SZ signal. Using a large chop throw to reduce the tendency to subtract SZ signal and chopping in azimuth so that sky rotation reduces the effect of point source contamination are both parts of this. Also, using the $`450`$$`\mu `$m array average to remove the spurious sky signal, thereby keeping SZ information on scales larger than the array size, maximizes the available signal in the data. Controls such as the blank fields, and to a lesser extent the reflector data, are required to check for systematic effects in both the instrument and analysis methods.

The second necessary component of this experiment is the use of custom-designed software to handle the data analysis, which is particularly important for fitting the differenced data directly to a model.

The third important component of this analysis is the removal of point sources in order to reduce their contaminating effect. Because sub-mm bright, high $`z`$ background sources are relatively common in cluster fields, observations of the SZ increment could suffer prohibitively from the effects of point sources. The basic principle we employ is that point sources affect only a few pixels and hence can be distinguished from an extended SZ effect. Our data have a confusion limit of approximately $`0.8`$mJy beam<sup>-1</sup> for SZ signals. Integration to SCUBA’s fundamental confusion limit of $`0.5`$mJy beam<sup>-1</sup> rms (which would require weeks of integration time) gives an SZ confusion of $`0.5`$mJy beam<sup>-1</sup>. Our simulations show negligible reduction in the confusion limit of these observations when point sources with amplitude less than about 5 mJy are removed (ZNH, ). This limit would require approximately 2 shifts of integration time. Unfortunately, decreasing the confusion limit floor below $`0.5`$ mJy beam<sup>-1</sup> would require higher angular resolution observations. A corollary to this is that SZ increment detections from other experiments with lower resolution than SCUBA, such as SuZIE and ACBAR, may be confusion dominated if unsupported by higher resolution data.

These measurements confirm the SZ effect in 2 galaxy clusters at an amplitude consistent with the measured decrements. A weak constraint on the kinetic effect in Cl$`\mathrm{\hspace{0.17em}0016}`$ is found which is consistent with the result in Benson et al. (2003). We quote a marginal 2$`\sigma `$ SZ increment detection in MS$`\mathrm{\hspace{0.17em}1054}`$, but cannot do much more with the data at this point.

It is clear that large reductions in the uncertainty of the value for the kinetic effect velocity in Cl$`\mathrm{\hspace{0.17em}0016}`$ will not be achievable because of source confusion with this type of instrument. This means that SCUBA is not capable of determining the amplitude of the kinetic effect in any specific cluster, given the expected level of SZ–derived peculiar velocities (e.g. Sheth & Diaferio 2001). It is likely that the only way to use SCUBA for constraining peculiar velocities is through a statistical survey of many clusters.

Nevertheless, we have demonstrated that, with sufficient care, it is possible to measure the $`850`$$`\mu `$m increment at the JCMT. Such SCUBA observations of the SZ effect at $`850`$$`\mu `$m can be combined with measurements at other frequencies to study individual clusters in more detail. Follow–up observations with high resolution sub-mm instruments are likely to become an important part of forthcoming SZ cluster surveys.

## Acknowledgments

We are grateful to Dr. Mark Birkinshaw for his helpful comments. This work was supported by the Natural Sciences and Engineering Research Council of Canada. The James Clerk Maxwell Telescope is operated by The Joint Astronomy Centre on behalf of the Particle Physics and Astronomy Research Council of the United Kingdom, the Netherlands Organization for Scientific Research, and the National Research Council of Canada. We would like to acknowledge the staff at JCMT for facilitating these observations. This research has made use of NASA’s Astrophysics Data System, the SIMBAD database, operated at CDS, Strasbourg, France, and the Canadian Astronomy Data Centre. E. P. was funded by NSERC and NASA grant NAG5-11489 during the course of this work.

This paper has been produced using the Royal Astronomical Society/Blackwell Science style file.


format: text

text:
# Relaxation kinetics in two-dimensional structures

## I Introduction

A large amount of work has been done in the domain of atomic structure kinetics since Burton, Cabrera and Frank presented, in a seminal paper, the first serious attempt to model the behavior of atoms adsorbed onto a vicinal surface bcf . Nevertheless, several problems remain unsolved. Many important advances in the techniques for the observation and manipulation of atoms now allow dealing with structures of decreasing size and, with the help of powerful computation, gaining insight into the basic mechanisms governing atomic dynamics. (For a review, see for example Ref. lag, ; zha, ; jen0,  and references therein).

In the last few years, considerable efforts have been dedicated to understanding how nanostructures relax. A classical description of structure shape equilibration was developed by Herring, Nichols and Mullins (HNM) mul . In this theory, the transport of mass responsible for relaxation occurs in a smooth and continous way via migration of border adatoms from regions of high curvature (high chemical potential $`\mu `$) to regions of low curvature (low $`\mu `$). It is assumed that structure border is rough enough to allow a continous (coarse-grained) description of it. However, as this is not valid for surfaces below the roughening temperature, the low temperature decay has been the subject of research both theoretically jen1 ; com ; tri  and experimentally ex1 ; ex2 .

Recently, studies of two-dimensional island shape relaxation with a simple model has revealed new, unexpected phenomena jen1 ; com . It was shown that two qualitatively different relaxation regimes exist. At high temperature, islands evolve toward equilibrium according to HNM. At low temperature, islands become faceted and a new driving mechanism appears. It was demonstrated that, in this case, shape relaxation occurs by steps, where the limiting process consists in the nucleation of new adatom rows on the flat edges of islands. In the above model adatoms lie on a triangular lattice and both equilibrium and kinetic properties depend on the single parameter $`ϵ_0/k_BT`$ (where $`k_B`$ is the Boltzmann constant and $`T`$ the temperature) since, for an adatom with $`n_i`$ nearest-neighbors (NN), $`n_iϵ_0`$ is both the potential energy and the kinetic barrier controlling migration.

This two-regime scenario is manifest in the dependence of the equilibration time $`t_{eq}`$ on temperature and island size $`N`$. At high temperature, $`t_{eq}N^2\mathrm{exp}(3\beta ϵ_0)`$ ($`\beta =1/k_BT`$), while at low $`T`$, $`t_{eq}N\mathrm{exp}(4\beta ϵ_0)`$. The regime to which relaxation of a given island belongs is determined by border roughness and depends not only on temperature but also on island size. In the $`TN`$ plane, the separation line between the two regimes is provided by the crossover island size $`N_c(T)\mathrm{exp}(\beta ϵ_0)`$, a rough indication of the size of the largest island that is fully faceted at $`T`$. Interestingly, if the temperature-dependent factors appearing in the leading terms of $`t_{eq}`$ are rewritten in terms of $`N_c`$, the properly-scaled equilibration time becomes a function of $`N/N_c`$ only, satisfying

$$\left(\frac{t_{eq}}{N_c^5}\right)\{\begin{array}{ccc}\left(\frac{N}{N_c}\right)& \text{for}& \frac{N}{N_c}1\\ & & \\ \left(\frac{N}{N_c}\right)^2& \text{for}& \frac{N}{N_c}1\end{array}\text{.}\text{}$$
(1)
A first problem of interest is the universality of this scaling law. The model described above is at best a first approximation to real situations and many details can be modified in order to improve it. For example, such a well-known effect as the Ehrlich-Schwoebel barrier schwo  can be accounted for by introducing new kinetic barriers in a phenomenological way; in an even more accurate approach, minimum energy transition paths could be calculated from interatomic forces neb . As a different set of kinetic barriers will entail, to the least, new activation energies, the question remains whether the above scaling behavior is a consequence of the model or, to the contrary, would carry over to more general situations.

Another closely-related problem is the relaxation of two-dimensional pores, i. e., islands of vacancies. In some sense, islands and pores are ‘specular images’ of one another. A given pore has the same boundary as the corresponding island, and so would be faceted or rough depending on whether the number of vacancies in it is smaller or larger than $`N_c(T)`$. It would be desirable to know if pore shape relaxation exhibits also a two-regime scenario and to analyze the possibility of scaling laws for the corresponding equilibration times.

As far as we know, the questions above have not been investigated, even in the case of simple models. In this work we study the main properties of shape relaxation in two-dimensional structures of adatoms. Our approach is two-fold: theoretical analysis and numerical simulations. In the analytic part, we treat the shape relaxation problem in a general way. Our study suggests the presence of two regimes, regardless of the specific set of kinetic barriers and, based on detailed-balance conditions, we show that similar scaling laws apply for both islands and pores. Since most of the arguments we invoke are heuristic, some independent confirmation of the derived properties is desirable. This is the goal of the numerical part. We analyze, using standard kinetic Monte Carlo (KMC) simulations, the relaxation of islands and pores according to two models considered in this paper, and compare the results with our theoretical predictions.

The paper is organized as follows. In Sec. II we discuss the causes for scaling behavior and derive the asymptotic laws of scaling functions. The models we use are defined in Sec. III and the outcome of the KMC simulations are presented in Sec. IV. Finally, we give our conclusions in Sec. V.

## II Analytic approach

In this section we investigate a generic model with adatoms lying on a triangular lattice where the total binding energy per NN pair is $`ϵ_0`$. Equilibrium properties depend on a single parameter, namely $`\beta ϵ_0`$, but kinetics involve in general a greater number of them. Our aim here is to establish scaling laws for both island and pore relaxation times.

Let us consider islands first. It is easy to see that many of the arguments presented in Ref. jen1, ; com,  remain applicable regardless of the specific set of kinetic barriers. First of all, $`N_c(T)`$, which separates rough and faceted islands at a given temperature, is the size at which the perimeter equals the average distance between border kinks. This does not depend on the transition barriers but only on the binding energies in the different configurations, so the form $`N_c\mathrm{exp}(\beta ϵ_0)`$ is still valid. Next, two qualitatively different modes of relaxation should exist. On the one hand, for islands greater that $`N_c`$, shape relaxation dynamics is adequately described by the HNM theory. It consists of a set of equations for the temporal evolution of the border curvature, where the kinetics is dominated by perimeter diffusion and always leads to a $`t_{eq}N^2`$ power law. We will call this regime rough relaxation mode (RRM). On the other hand, mechanisms driving islands smaller than $`N_c`$ toward equilibrium are also quite general. In this scenario, islands spend most of the time in highly-faceted configurations. We will call this regime faceted relaxation mode (FRM). Evolution occurs by steps consisting in nucleation and stabilization of new adatom rows. The time associated to each of these steps is proportional to facet length ($`N^{1/2}`$) and, given an initial configuration, a number $`N^{1/2}`$ of rows needs to be created in order to attain the equilibrium shape; thus $`t_{eq}N`$.

The size exponents for equilibration time thus appear to be universal (1 in the FRM, 2 in the RRM). The activation energies however do depend in general on the kinetic barriers. The activation energy in the FRM corresponds to the energy characteristic of row nucleation, while in the RRM it will be equivalent to the characteristic energy for diffusion along the island border. If we call these energies $`E_F`$ and $`E_R`$ for the FRM and RRM, respectively, the arguments presented so far allow us to write the ansatz

$$t_{eq}\{\begin{array}{ccc}\left(\frac{N}{N_c}\right)e^{\beta (E_F+ϵ_0)}& \text{for}& \frac{N}{N_c}1\\ & & \\ \left(\frac{N}{N_c}\right)^2e^{\beta (E_R+2ϵ_0)}& \text{for}& \frac{N}{N_c}1\end{array}\text{.}\text{}$$
(2)
At this point, it is clear that the occurrence of scaling is limited to the cases where

$$E_F=E_R+ϵ_0\text{;}$$
(3)
otherwise, the exponential factors in Eq. (2) would be different and no natural time scale \[such as $`N^5`$ in Eq. (1)\], allowing scaled $`t_{eq}`$ to be expressed as a function of $`N/N_c`$, would exist. Let us, then, analyze the connexion between the two energies.

As noted above, the limiting step in the RRM is adatom border diffusion. This is in turn limited by the elementary processes sketched in Fig. 1, in which adatoms pass from one kink to another by jumping around a corner; thus we can readily equal $`E_R`$ and the energy of this process.

For the FRM, the dynamics is determined by row nucleation, which occurs when two adatoms meet on a flat island edge. The nucleation rate may thus be calculated by the product of two factors: the rate of adatoms entering a facet and the probability that another adatom lies on the same facet. In this regime, the main sources of adatoms are facet ends, where the most weakly bound adatoms are found. Those adatoms move onto flat facets by basically the process sketched in Fig. 1 (A); the first factor is thus $`\mathrm{exp}(\beta E_R)`$. The second factor is $`\mathrm{exp}(\beta ϵ_0)`$ because the system gains an energy $`ϵ_0`$ when an adatom moves from a kink to a flat edge. Thus, we have $`E_F=E_R+ϵ_0`$, i. e., Eq. (3), and the scaling form

$$\frac{t_{eq}}{N_c^\alpha }\{\begin{array}{ccc}\left(\frac{N}{N_c}\right)& \text{for}& \frac{N}{N_c}1\\ & & \\ \left(\frac{N}{N_c}\right)^2& \text{for}& \frac{N}{N_c}1\end{array}\text{,}\text{}$$
(4)
where the exponent $`\alpha `$ is

$$\alpha =2+\frac{E_R}{ϵ_0}\text{.}\text{}$$
(5)
We consider now the case of pores. Since, as mentioned before, a pore will also be rough or faceted depending on whether it is larger or smaller than $`N_c`$, it is worth analyzing each situation separately.

For the relaxation of rough pores, the solution is immediate: the kinetic path is the same as that followed by the corresponding island (via a vacancy-atom exchange mechanism). This pore-island symmetry is clearly evident in the HNM theory (see for example Ref. jen1, ), which involves only border curvature and perimeter atomic diffusion and hence does not depend on the kind of particles enclosed within the border. As a corollary, the equilibration time for pores should also satisfy the asymptotic law established by the second equation in Eq. (2); this will be verified numerically in Sec. IV.

For faceted pores, it is clear that relaxation will occur, as in the case of faceted islands, by steps consisting of the nucleation and stabilization of new rows. However, despite this qualitative analogy, the two processes are not equivalent but, rather, complementary; adatom migration occurs from short to long facets in islands, while the opposite is true in pores. For the latter, row nucleation occurs when, as shown schematically in Fig. 2, an adatom that lies on a long facet moves to an internal corner of the pore. The new row will grow with the arrival of more adatoms from the now ‘open’ long facet, and become stable when fully filled.

The similarity of the relaxation mechanisms allows us to apply, in the case of pores, the same kind of arguments we invoked for the dependence of $`t_{eq}`$ on the size of faceted islands, which again yields the power law $`t_{eq}N`$. In order to calculate the activation energy –which corresponds to the nucleation process illustrated in Fig. 2– we consider first the inverse process, i. e., the migration of an adatom from site $`\mathrm{𝟐}`$ to site $`\mathrm{𝟏}`$. In this case, the adatom first leaves a kink, then diffuses along a flat border and, finally, goes into an adjacent row by jumping around a corner. This is similar to the process sketched in Fig. 1 (B) and must therefore have an activation energy $`E_R`$. Knowledge of this last energy allows us to find that correponding to nucleation. They are related by the detailed-balance condition which, taking into accout that the binding energy difference between sites $`\mathrm{𝟏}`$ and $`\mathrm{𝟐}`$ is $`ϵ_0`$, finally leads to the result that the activation energies for pore shape relaxation are also connected by relation (3).

As an interesting consequence, the equilibration time for islands or pores, when properly scaled with $`N_c^\alpha `$, becomes a function of $`N/N_c`$ only that satisfies the asymptotic rules (4). We stress that even though the HNM theory implies an equivalent relaxation mode for both kinds of structures in the RRM, the same symmetry cannot be expected to hold in the FRM. The asymptotic forms give the exponents of the leading terms but different prefactors may appear. In general, there will be one scaling function corresponding to islands and another scaling function for pores. For very small values of their arguments ($`NN_c`$) they will be parallel (on a log scale) but not necessarily equal, while at the other extreme ($`NN_c`$) the curves should collapse as required by HNM theory. Using KMC we will see, in Sec. IV, that these properties of the scaling functions are verified.

## III The Models

In order to asses the analytical predictions above, we investigated two different models, denoted I and II, which have the following common characteristics: The substrate is represented by a triangular lattice and the total binding energy per NN pair is $`ϵ_0`$. Adatom hoppings are only possible between NN sites and a given adatom in an initial site $`a`$ jumps to a final empty site $`b`$ with a probability per unit time $`P_{ab}=\nu \mathrm{exp}(\beta E_{ab})`$, where $`\nu `$ is a constant frequency and $`E_{ab}`$ the kinetic barrier. In order to avoid detachment, the jumps are forbidden if the number of NN adatoms in the final site is 0. As an additional simplification, the motion of adatoms with initially more than 4 NN is also forbidden. This approximation is justified because of the high energy barriers of the corresponding processes; it is introduced in order to accelerate the simulations. Differences between models I and II lie only in the set of kinetic barriers.

For Model I, which is the same as the one used in Ref. jen1, ; com, , $`E_{ab}`$ depends only on the number $`n_a`$ of NN adatoms in the initial site $`a`$. It takes the value $`E_{ab}=n_aϵ_0`$ when $`1<n_a<5`$ and $`E_{ab}=ϵ/10`$ when $`n_a=1`$, regardless of the number of NN adatoms in $`b`$.

For Model II, we use a more accurate set of activation energies, obtained from a detailed study of the kinetic barriers for adatoms lying on a triangular lattice and interacting via a Lennard-Jones potential thibault . It is useful to define this model with the help of Fig. 3, where the different barriers are represented. In the top part of the figure, adatoms are represented by circles while crosses indicate possible binding sites. The bottom part shows the corresponding energy landscape, i. e., kinetic barriers separating adjacent local minima. Thus, in this model, we take into account the final-state configuration when evaluating the transition-state energy. (Note however that the height of the barrier depends on the number of NN in the initial and final states but not on their precise position). Model II thus incorporates such important features as the Ehrlich-Schwoebel barrier. This is apparent if we look at Fig. 3: the barrier from $`\mathrm{𝟐}`$ to $`\mathrm{𝟏}`$ is lower than the barrier from $`\mathrm{𝟐}`$ to $`\mathrm{𝟑}`$; the barrier from $`\mathrm{𝟒}`$ to $`\mathrm{𝟑}`$ is higher than the barrier from $`\mathrm{𝟒}`$ to $`\mathrm{𝟓}`$, etc. The final parameters in this model are $`E_1=0`$, $`E_2=0.5ϵ_0`$, $`E_3=ϵ_0`$, $`E_4=1.5ϵ_0`$, $`E_5=2ϵ_0`$, $`E_6=4ϵ_0`$. In order to facilitate comparison, Table 1 lists, for both models, the kinetic barriers corresponding to a given elementary jump as a function of the number of NN before ($`n_i`$) and after ($`n_f`$) the hopping.

Before closing this section, we calculate the activation energies for each model as derived from the analysis of Sec. II. Remember that $`E_R`$ is the energy characteristic of the process sketched in Fig. 1 (A). For Model I, this process corresponds to an adatom with 3 NN before the jump, so $`E_R^I=3ϵ_0`$, and, using Eqs. (2) and (5), $`E_F^I=4ϵ_0`$ and $`\alpha ^I=5`$, as obtained in Ref. jen1, ; com,  (but not yet tested with pores). For Model II, the process in question is (for example) the transition from $`\mathrm{𝟒}`$ to $`\mathrm{𝟐}`$ in Fig. 3, which leads to $`E_R^{II}=2ϵ_0`$, $`E_F^{II}=3ϵ_0`$ and $`\alpha ^{II}=4`$.

## IV Numerical results

As mentioned in Sec. I, a numerical confirmation of the results of Sec. II, obtained in a heuristic manner, is in order. In this section, we do this using KMC simulations to explore shape relaxation according to models I and II. We consider islands and pores with sizes between 400 and 6500 particles and perform, for each model, standard KMC simulations. The values of $`\beta `$ (in $`1/ϵ_0`$ units) range between 2 and 17. The results shown in this section correspond to averages over a number of samples, varying from 4 for the largest structures to 30 for the smallest ones.

Following Ref. jen1, ; com, , we use the aspect ratio $`\alpha `$, defined as the ratio of the $`x`$ and $`y`$ gyration radii, to characterize the state of the structures. We start each simulation with an aspect ratio around 10 and define the equilibration time $`t_{eq}`$ as the time at which $`\alpha `$ first becomes less than 1. First of all, we analyze the dependence of $`t_{eq}`$ on $`N`$ (for fixed $`\beta `$) and, separately, on $`\beta `$ (for fixed $`N`$). Later in this section, we will check the scaling laws.

Let us start with Model I. In order to test the behavior of $`t_{eq}`$ as a function of $`N`$, we plot in Fig. 4 ($`a`$) $`\mathrm{ln}t_{eq}`$ against $`\mathrm{ln}N`$ for several temperatures. Filled and empty symbols correspond to islands and pores, respectively. At high (low) enough temperature, $`t_{eq}`$ is expected to scale as $`N^2`$ ($`N`$). This scaling behavior is represented in Fig. 4 ($`a`$) by the lines with slopes of 2 (lower) and 1 (upper), which bracket the KMC results corresponding to the highest and lowest temperatures used in our simulations, respectively. Using these lines as a guide, it is indeed clear that the the size exponent changes from 1 to 2 when the temperature is increased, consistent with the numerical findings. Note that the agreement is as good for pores as it is for islands.

We check now our predictions for the activation energies. To do so, in Fig. 4 ($`b`$), $`\mathrm{ln}t_{eq}`$ is plotted against $`\beta `$ for different sizes. The asymptotic behaviors are represented by the straight lines which have slopes of $`3ϵ_0`$ (low $`\beta `$) and $`4ϵ_0`$ (high $`\beta `$). The KMC data again show excellent agreement with the predicted behavior, to within a trivial prefactor (i. e., slopes agree). Let us remark that the results for island relaxation are given here to allow comparison with the (new) results for pores. We note however a small discrepancy between our numerical results and those of Ref. com, : although the size exponents and activation energies coincide, our equilibration times are found to be lower (by a factor of $`10`$) than those reported in Ref. com, .

The corresponding results for Model II are presented in Fig. 5. Note that, also in this case, the theory-simulation agreement is good for both types of structures. Thus, Figs. 4 and 5 confirm our predictions that, on the one hand, the size exponents are universal (1 and 2) and, on the other hand, $`E_R`$ and $`E_F`$ depend on the specific set of microscopic kinetic barriers.

Finally, in order to test the specific scaling laws, we plot, on a log-log scale, $`(t_{eq}/Nc^5)`$ against $`(N/Nc)`$ for Model I (Fig. 6) and $`(t_{eq}/Nc^4)`$ against $`(N/Nc)`$ for Model II (Fig. 7). We have used $`\alpha ^I`$ and $`\alpha ^{II}`$ as calculated at the end of Sec. III. According to the analysis presented in Sec. II, $`t_{eq}/N_c^\alpha `$ should scale with exponents 1 and 2 for very small and very large values of its argument ($`N/N_c`$). In Figs. 6 and 7, we have drawn lines with slopes of 1 and 2 to indicate these asymptotic behaviors, which are evidently closely verified. For concreteness, in these plots, we use the same $`N_c=.25\mathrm{exp}(\beta ϵ_0)`$ as in Ref. jen1, ; the exact value of the geometrical prefactor should not affect the final result in a significant manner.

It is interesting to note that, besides confirming the expected asymptotic power laws, in both figures all the points corresponding to islands collapse on one curve, and all the points corresponding to pores collapse on another one (depending via $`\alpha `$ on the microscopic details of the model), thus giving additional support to the scaling ansatz. Furthermore, although for each model the island and pore scaling functions are different in the FRM, they become equivalent in the RRM, as required by the perfect pore-island symmetry of relaxation in the latter regime.

## V Conclusions

By considering a generic model of adatoms lying on a triangular lattice, we have studied the problem of shape relaxation of islands and pores. The arguments employed in this work allows the main properties of the equilibration time $`t_{eq}`$ to be calculated as a function of temperature and size. Because our arguments are somewhat heuristic, KMC simulations, using two different kinetic models, were also carried out. The numerical results confirm our theoretical predictions. For both islands and pores two qualitatively different modes of relaxation (FRM and RRM) are found, as well as the line $`N_c(T)`$ that separates these regimes in the temperature-size plane. It was shown that, although size exponents are universal (1 in FRM, 2 in RRM), the activation energies corresponding to each mode depend on the microscopic details of the kinetic model. Scaling behaviors were found for the equilibration time: the properly scaled $`t_{eq}`$ depends also, via $`\alpha `$, on the set of elementary kinetic barriers of the model. Furthermore, for a given model, the specific scaling function for islands is in general different from the one corresponding to pores. While in the FRM both functions scale with the same exponent ($`=1`$) as a consequence of the detailed-balance condition, in the RRM the two functions collapse because of the perfect island-pore symmetry of HNM theory mul .

## AKNOWLEDGMENTS

We are grateful to N. Combe and P. Jensen for interesting discussions. We thank F. Bussières, M. de la Chevrotière, J. Richer, M. Simard and C. Hudon for help with the numerical calculations, and P. Thibault who provided the data needed for Model II. This work was supported by grants from the Natural Sciences and Engineering Research Council (NSERC) of Canada and the “Fonds Québécois de la recherche sur la nature et les technologies” (FQRNT) of the Province of Québec. We are indebted to the “Réseau québécois de calcul de haute performance” (RQCHP) for generous allocations of computer resources.


text:
# An accretion disk model for periodic timing variations of pulsars

## 1 Introduction

Two isolated pulsars (PSR B1828$``$11 & B1642$``$03) present the phenomenon of long$``$term, highly periodic and correlated variations in both the pulse shape and the slow$``$down rate (for PSR B1828$``$11, its variations are best described as harmonically related sinusoids, with periods of approximately 1,000, 500 and 250 days, see Stairs et al. 2000; for PSR B1642$``$03, its timing residuals exhibit cyclical changes with amplitude varying from 15 to 80 ms and spacing of maxima varying from 3 to 7 yr, see Cordes 1993 and Shabanova et al. 2001). This phenomenon provides the most compelling evidence for precession (Link & Epstein 2001). Some authors (Jones & Andersson 2001, Link & Epstein 2001, Rezania 2002) proposed different models to explain this phenomenon within the framework of free precession. However, further investigations are expected on whether the force of vortices pinning to the stellar crust is strong enough to damp out the free precession of neutron stars (see the conclusion section for details).

The pulsating X$``$ray binary Her X$``$1 also shows an unusual kind of periodic timing behavior in the form of an on$``$off cycle about 35 days. Sarazin et al. (1980) showed that a moderately massive disk ($`10^510^2`$ $`M_{}`$) can cause the precession of the central compact star. Trumper et al. (1983) used the free precession of the neutron star to explain the 35$``$day period. Qiao & Cheng (1989) suggested that instead of the free precession of the neutron star, the periodic variation is due to the precession of the accretion disk driven by a normal star.

In this paper, we present an accretion disk model to account for the precession of pulsars. The total torque exerted on the pulsar by the fossil disk, which causes the pulsar to precess, is nonzero, due to the asymmetry in pulsar shape and nonalignment between the disk rotation axis and the axis of pulsar moment of inertia. Hence, we present this forced precession model (disk$``$induced precession), other than the free precession model.

In fact, such a system, i.e. a pulsar with an accretion disk, has been invoked to describe the nature of two kinds of isolated X$``$ray pulsars: the anomalous X$``$ray pulsars (AXPs, which have persistent X$``$ray pulsations but few burst) and the soft gamma$``$ray repeaters (SGRs, which have multiple bursts of soft gamma$``$rays) (see Xue et al. 2003 for a review). The most remarkable difference between these two kinds of radio$``$quiet X$``$ray pulsars and the radio pulsars is that their X$``$ray luminosities (usually $`L_\mathrm{X}`$ $``$ 1035 ergs/s) are several orders of magnitude higher than their rotational energy loss rates. To explain this and other exotic facts, two alternative models, viz. the magnetar models (in which the pulsar is suggested to be a neutron star with super$``$strong magnetic field) and the accretion disk models, were proposed (see Xue et al. 2003 for a review). Under the framework of the accretion disk model, a problem would come if the precession of radio pulsars (e.g. PSR B1811$``$11 and B1642$``$03) was really caused by the accretion disk: why do some pulsars with disks behave as ordinary radio pulsars while others behave as AXPs/SGRs? In this paper, we discuss this problem in section 3 after presenting the detailed calculation in section 2.

## 2 Pulsar precession caused by an accretion disk

### 2.1 Geometry and parameters

Consider the neutron star as a rotation ellipsoid with the principal moment of inertia $`I_1=I_2I_3`$ and the corresponding radii $`a=ab`$. The geometry of the system is shown in Fig.1: $`O`$$``$$`xyz`$ is the coordinate system fixed with the neutron star, $`Ox`$, $`Oy`$ and $`Oz`$ are the principal axes of moment of inertia of the neutron star, the axis $`OA`$ is the rotation axis of the neutron star (for simplicity, we assume the axis $`Oz`$ is aligned with the axis $`OA`$ in this paper; in this case, free precession doesn’t exist and we can see how the disk makes the neutron star precess more clearly), $`O`$$``$$`\xi \eta \zeta `$ is the coordinate system fixed with the disk and the disk plane is $`O\xi \eta `$, $`\theta `$, $`\phi `$ and $`\varphi `$ are Euler angles, $`\omega `$ is the angular velocity of neutron star rotation. We will derive the formulae of the gravitational potential in the neutron star$``$disk system and the angular velocity of neutron star precession. And then, as an example, we model the precession period of PSR B1818$``$11. The parameters about PSR B1828$``$11 are as follows: the rotational period of this pulsar $`P=0.405`$s ($`\omega =2\pi /P`$), the period derivative $`\dot{P}=6.0\times 10^{14}`$ ss-1, the precession period $`P_{\mathrm{pre}}1000`$ days, the magnetic field $`B=5.0\times 10^{12}`$G (Stairs et al. 2000).

### 2.2 Equations and numerical results

The gravitational potential between the neutron star and the disk is

$`V`$ $`=`$ $`{\displaystyle \frac{2GMM_0}{c+d}}`$ (1)
$`\times `$ $`[1{\displaystyle \frac{b^2+2a^2}{10cd}}+{\displaystyle \frac{3b^2}{20cd}}\mathrm{sin}^2\theta +{\displaystyle \frac{3a^2}{20cd}}(1+\mathrm{cos}^2\theta )],`$
where $`M`$, $`a`$ and $`b`$ are the mass, the equatorial and longitudinal radii of the neutron star, respectively, $`M_0`$ is the mass of the disk, $`c`$ and $`d`$ are the inner and outer radii of the disk, respectively.

The moment of inertia of the neutron star is $`I_\mathrm{z}=(8/15)\pi \rho a^4b`$, where $`\rho =3M/(4\pi a^2b)`$ is the density of the neutron star. Since the rotation angular velocity of the neutron star $`\omega `$ is much greater than the precession angular velocity of the neutron star $`\dot{\varphi }`$, the angular momentum of the neutron star can be written as $`L_0I_\mathrm{z}\omega `$. According to the Lagrange’s equation $`L_0\dot{\varphi }\mathrm{sin}\theta =V/\theta `$ (Qiao & Cheng 1989), we get the precession angular velocity of the neutron star:

$$\dot{\varphi }=\frac{3GM_0\mathrm{cos}\theta }{2cd(d+c)\omega }(1\frac{b^2}{a^2}).$$
(2)
In terms of the study on the stability of the Maclaurin spheroids, we can obtain the oblateness of this pulsar:

$$ϵ=\frac{ab}{a}=2.27\times 10^6.$$
(3)
For simplification, we assume $`cdR`$, so that

$$\dot{\varphi }=\frac{3GM_0\mathrm{cos}\theta }{4R^3\omega }[1(1ϵ)^2].$$
(4)
Fig.2 shows the relations between $`M_0`$ and $`R`$ derived from the 1000$``$day precession period for a group of $`\theta `$. In principle, $`M_0`$ may vary from less than $`0.001M_{}`$ to even $`0.1M_{}`$ (The amount of the fall$``$back material ($`M_{\mathrm{fb}}`$) which forms the disk after supernova explosion is uncertain, however, $`M_{\mathrm{fb}}0.1M_{}`$ is supposed to be reasonable, see Chatterjee et al. 2000), and $`R`$ from about $`1R_{\mathrm{co}}`$ (or smaller) to $`23R_{\mathrm{co}}`$ when $`\theta `$ varies within the range of $`090^o`$, where $`R_{\mathrm{co}}`$ is the corotation radius of the disk with the pulsar. Thus, the precession period of the PSR B1828$``$11 is well elucidated.

## 3 Possible evolutionary scenario in the pulsar$``$disk system

One may ask where the accretion disk comes from. Usually, a disk around the pulsar can be formed: 1. by accreting the matter from a companion star in the binary system (approach I); 2. by the fall$``$back of the material after a supernova explosion (approach II). In the former process, there are two ways of mass transfer: matter flowing through the first Lagrange point due to the critical Roche Lobe abounding with the companion matter, or by the stellar wind (Nagase 1989).

An interesting issue is that what the pulsar in the pulsar$``$disk system will behave as, radio pulsar with precession or AXP/SGR? In the accretion model of X$``$ray binary, accreting process may experience some “phases”, e.g. “radio pulsar phase”, “sleeping phase” , “propeller phase” and “accretion phase” (Illarionov & Sunyaev 1975, Henriches 1983). The case is suggested to be similar for a pulsar$``$disk system formed through either approach I or II. The PSR B1828$``$11, B1642$``$03, AXPs and SGRs are all isolated “neutron stars”. What are the links between them? Here we present a possible series of evolutionary tracks for the pulsar$``$disk system (cf. Chatterjee et al. 2000 and Henriches 1983).

(a) After the supernova explosion, the newborn neutron star rotates so fast that its magnetospheric radius $`R_\mathrm{m}`$ is much larger than its light cylinder radius $`R_{\mathrm{LC}}`$. The coupling between the disk and the neutron star is unlikely. The young pulsar initially emits a very strong relativistic wind together with highly energetic magnetic dipole radiation. The pulsar wind and photons blow a cavity in the fall$``$back flow and no matter can accrete, as is similar in the binary system (Henriches 1983). Thus the neutron star may live in the “radio pulsar phase”, i.e. the radio emission can propagate away.

(b) The neutron star slows down as its rotational energy reduces. Its radiation and the high energy particle flow become weakened, thus the cavity decreases and the radio emission is finally quenched. When the “radio pulsar phase” ends, accretion down to the surface of the neutron star is not expected immediately. The so$``$called “sleeping phase” takes place. This phase lasts until the pulsar rotation slows down enough to make accretion possible.

(c) Subsequently, the neutron star enters the “propeller phase” since accretion onto the neutron star is inefficient due to centrifugal forces acting on the matter. As the rotation slows down further, accretion will be permitted and the neutron star will behave as an accretion$``$driven X$``$ray source. Gradually, the system enters a quasi$``$equilibrium “tracking phase” (happening at the time of $`t_{\mathrm{trans}}`$), in which the spin of the star roughly matches the rotation of the disk at $`R_\mathrm{m}`$. In this phase, the neutron star will be X$``$ray bright due to accretion. It is possible that the advection$``$dominated accretion flow (ADAF) will occur at the time of $`t_{\mathrm{ADAF}}`$, when much of the mass will be ejected prior to reaching the neutron star surface. The X$``$ray luminosity of the neutron star will decline rapidly with time due to $`\dot{M}_\mathrm{X}`$ (the rate of accretion onto the star surface) diminishing. Therefore, the bright AXPs may exist only between $`t_{\mathrm{trans}}`$ and 2$`t_{\mathrm{ADAF}}`$ (about 5000 $``$ 40000 years).

(d) Finally, accretion ceases since the matter of the disk is exhausted. Thus no accretion$``$driven X$``$ray pulsar can be observed and neither can radio emission be emitted owing to the long rotational period of the neutron star (Qiao & Zhang 1996, Zhang et al. 2000). Consequently, the neutron star may be observed as a dim isolated neutron star, a weak soft X$``$ray source powered by rotational energy loss or thermal radiation.

## 4 Conclusions and discussions

The accretion disk model for periodic timing variations of pulsars (e.g. PSR B1828$``$11) can elucidate the observed precession period. According to the analysis of the possible evolutionary scenario discussed above, the PSR B1828$``$11 may be just in the “radio pulsar phase” before the “sleeping phase” and the “propeller phase”. If our model proves to be true, the link between radio pulsars and AXPs may be established and the fossil disk model for AXPs may be strengthened. It also suggests that the neutron star$``$fossil disk system might exist. Some discussions are listed below.

1. Sarazin et al. (1980) pointed out when the angular momentum of the disk is larger than that of the neutron star, the neutron star will precess due to the effect of the accretion disk and the disk will maintain its orientation. This is just our case for the precession of the PSR B1828$``$11.

2. Why is the radio emission expected in our model but not expected in the case of close X$``$ray binary system? Illarionov & Sunyaev (1975) suggested that the large optical depth of the stellar wind for free$``$free absorption hinders the detection of radio pulsars. In our model, a fossil disk is formed and there is no enough stellar wind to absorb the radio emission. A case in point is the PSR B1259$``$63, which has a $`B_\mathrm{e}`$ companion star with high ellipticity and long orbital period. This neutron star is observed as a radio pulsar when the companion encircles far away from it. But no radio emission is detected when the companion approaches very near the pulsar (Johnston et al. 1992).

3. The condition when the pulsar radio emission quenches is not clear yet. A possible case may be that the cavity radius becomes smaller than the light cylinder radius (Henriches 1983). In our model, although the corotation radius of the PSR B1828$``$11 $`R_{\mathrm{co}}=(GM/\omega ^2)^{1/3}9.19\times 10^7`$cm is one order of magnitude smaller than its light cylinder radius $`R_{\mathrm{LC}}=c/\omega =1.93\times 10^9`$cm, and the disk radius $`R`$ approaches $`23R_{\mathrm{co}}`$, this neutron star is still in the “radio pulsar phase”. This presents an example for further investigations.

4. We have ignored the angular momentum of the pinned superfluid in the dynamic equations of the calculation. Can the pinning force be strong enough to damp any precession of neutron stars? In principle, detailed computations may determine whether the vortices are pinned, but a definitive conclusion on the nature has not been reached yet due to various uncertainties in the microscopic physics. Based on the calculation of Shaham (1977) and Sedrakian et al. (1999), Stairs et al. (2000) conclude that the torque$``$free precession of PSR B1828$``$11 should be damped out due to the dissipation caused by the vortices pinning to the stellar crust on the timescale of several hundred precession periods. Nevertheless, a further consideration of Shaham (1986) indicates an unpinned superfluid in the case of precession. It is also noted by Link & Cutler (2002), that hydrodynamic forces present in a precessing star are probably sufficient to unpin all of the vortices of the inner crust. In addition, the 35$``$d periodicity of the accreting binary system (Her X$``$1) has been interpreted as precession (Shakura, Postnov & Prokhorov 1998, and the references therein), which could be another example of forced precessions. We have therefore a tendency to suggest that, the observation of PSR B1828$``$11 precession may imply that the vertex pinning in this pulsar is much weaker than that predicted previously, at least it is weaker than the possible disk$``$driven precession effect.

5. In fact, the phenomenon of the secular periodic timing variations, which also exists in many other radio pulsars (e.g. Crab and Vela pulsar), has been investigated both observationally and theoretically (Sedrakian et al. 1999 and the references therein). If a neutron star were a rigid or semirigid solid star, it would precess (Pines & Shaham 1974); as is mentioned above, the interaction between the superfluid vortices and the stellar crust needs further considerations. Ruderman (2001) pointed out that this is a crucial problem which will raise doubts about canonical descriptions of neutron star structure until it is solved. No matter the precession is torque$``$free or forced, the neutron star structure still needs to be investigated more. Instead of neutron star, Xu (2003) proposed a solid strange star model, which may be possible and also needs to be investigated more.

6. Our model cannot exclude the possibility of a companion star (instead of a fossil disk) which can produce the precession phenomenon as well. In this case, two limits should be considered: the Roche Limit $`d_\mathrm{R}`$ (if the companion approached the neutron star closer than this limit, it would be torn asunder by tidal forces) and the Instability Limit $`d_\mathrm{I}`$ (if the company strayed farther than this limit, it would escape due to the differential perturbations from other celestial bodies). The details of this case will be studied in another paper.

###### Acknowledgements.

We are grateful to K. J. Lee et al. for their helpful discussions. This work is supported by National Nature Science Foundation of China (10173002, 10073001).


text:
# Anomalous heat conduction and anomalous diffusion in one dimensional systems

## Abstract

We establish a connection between anomalous heat conduction and anomalous diffusion in one dimensional systems. It is shown that if the mean square of the displacement of the particle is $`\mathrm{\Delta }x^2=2Dt^\alpha (0<\alpha 2)`$, then the thermal conductivity can be expressed in terms of the system size $`L`$ as $`\kappa =cL^\beta `$ with $`\beta =22/\alpha `$. This result predicts that a normal diffusion ($`\alpha =1`$) implies a normal heat conduction obeying the Fourier law ($`\beta =0`$), a superdiffusion ($`\alpha >1`$) implies an anomalous heat conduction with a divergent thermal conductivity ($`\beta >0`$), and more interestingly, a subdiffusion ($`\alpha <1`$) implies an anomalous heat conduction with a convergent thermal conductivity ($`\beta <0`$), consequently, the system is a thermal insulator in the thermodynamic limit. Existing numerical data support our results.

Does heat conduction in one dimensional (1D) systems obey the Fourier law? If it does, what are the necessary and sufficient conditions? If it does not, then what is the reason and how the thermal conductivity diverges/converges with the system size $`L`$? These questions have attracted increasing attention in recent years . Although some progress have been achieved, many puzzles remain. For example, in an attempt to establish a connection between heat conduction and the underlying microscopic dynamics, there exist some controversial examples. In the ding-a-ling model Casati et al show that as onset of global chaos the heat conduction crosses over from an abnormal one to a normal one obeying the Fourier law. It is thus concluded that the chaos is a deciding factor. Later on, in order to show the exponential instability is a necessary condition, Alonso et al  studied the heat conduction in a Lorentz gas channel, a quasi 1D billiard with circular scatterers, and found that the heat conduction obeys the Fourier law. However, the results from 1D Ehrenfest gas channels , in which the Lyapunov exponent is zero, show that the Fourier heat law might not have any direct connection to the underlying dynamical chaos, because the heat conduction can be normal and abnormal, depending on whether or not the disorder is introduced.

Recently, a quasi 1D triangle billiard model, which consists of two parallel lines of length $`L`$ at distance $`d`$ and a series of triangular scatterers, has been introduced and studied. In this model, no particle can move between the two reservoirs without suffering elastic collisions with the triangles. Therefore this model is analogous to the Lorentz gas channel studied in  with triangles instead of discs and the essential difference is that in the triangular model the dynamical instability is linear and therefore the Lyapunov exponent is zero. It is found that the motion inside the irrational triangle channel (the internal angles are irrational multiples of $`\pi `$) is diffusive and has a normal heat conduction. Therefore deterministic diffusion and normal heat transport which are usually associated to full hyperbolicity, can take place in systems without exponential instability. Another example is the Fermi-Pasta-Ulam (FPU) model which has non-zero Lyapunov exponent, however the heat conduction in this model does not obey the Fourier law.

The heat conduction in the rational triangle model (the internal angles are rational multiples of $`\pi `$) and in the FPU model is anomalous and does not obey the Fourier law, the thermal conductivity $`\kappa `$ diverges with system size $`L`$ as $`L^\beta `$ with $`\beta =0.22`$ for the rational triangle model, and $`0.34<\beta <0.44`$ for the FPU model. Indeed, similar divergent behavior has been observed in many 1D systems. For example, in the binary hard sphere model, $`0.22<\beta <0.35`$, in single wall nanotubes $`0.22<\beta <0.37`$ , and in many classical lattices such as the harmonic lattice, $`\beta =1`$, disordered harmonic lattice, $`\beta =1/2`$ , and the Frenkel-Kontorova (FK) model under the condition of $`T/K>>1`$, $`0<\beta <1`$, where $`T`$ is temperature, and $`K`$ is the effective amplitude of a sinusoidal on-site potential.

Obviously, a universal value of $`\beta `$ does not exist, it differs from model to model. Most recently, Narayan and Ramaswamy show theoretically that in a 1D momentum-conserving continuous system, the heat conduction is anomalous, and the thermal conductivity diverges with system size $`L`$ as $`L^{1/3}`$. Up to now, in all available numerical results only the heat conduction in a (5,5) single wall nanotubes shows an exponent ($`\beta 0.32`$) close to this 1/3. Despite the conduction mechanism is similar, a (10,10) single wall nanotube shows different value for unknown reasons. The numerical results from other models such as the FPU model, the harmonic model and other billiards models deviate largely from this value for reasons to be investigated.

On the other hand, even if the momentum conservation breaks down, the heat conduction can be anomalous such as that in the Frenkel-Kontorova model. The question becomes how to explain this anomalous heat conduction, in particular the value of the exponent $`\beta `$ in the thermal conductivity. A general theory is still lacking. The only existing theory is for the 1D harmonic chain, in which the phonons transport along the chain ballistically and the thermal conductivity, $`\kappa `$, diverges as $`L`$, i.e. $`\beta =1`$.

In this Letter, we would like to find a microscopic origin of the anomalous heat conduction observed in many 1D models. We shall not restrict to any specific model. This should give us a more general way to understand the heat conduction in 1D systems.

As is well known that, depending on the value of exponent $`\alpha `$ in the mean square of displacement of the particle, $`\mathrm{\Delta }x^2=2Dt^\alpha `$ with $`0<\alpha 2`$, 1D microscopic motion can be classified into ballistic motion, $`\alpha =2`$, superdiffusion, $`1<\alpha <2`$, normal diffusion, $`\alpha =1`$, and subdiffusion $`\alpha <1`$. Ballistic transport is observed in the harmonic lattice. Normal diffusion shows up in the FK model in a certain parameter regime, the disordered FPU model, the Lorentz gas channel, the disordered Ehrenfest gas channel, the irrational triangle channel, and the alternative mass hard-core potential model. In some billiard models, superdiffusion and subdiffusion are observed. Superdiffusion and subdiffusion can be studied from the fractional Fokker-Planck equation, for detailed theoretical investigation and discussion about the anomalous diffusion, please refer to review articles and the references therein.

To establish a connection between the microscopic process and the macroscopic heat conduction, let’s consider a 1D model of length $`L`$ whose two ends are put into contact with thermal baths of temperature $`T_L`$ and $`T_R`$ for the left end and the right end, respectively. Suppose the energy is transported by energy carriers (they are phonons in lattices and particles in billiard channels) from left heat bath to the right heat bath and vice versa. If the mean square of displacement of the carrier, with velocity $`v`$, inside the system can be described by,

$$\mathrm{\Delta }x^2=2Dv^\alpha t^\alpha $$
(1)
then the so-called “mean first passage time” (MFPT) is

$$t_{LR}=\frac{4\gamma }{\alpha \pi v}\left(\frac{2L}{\pi \sqrt{D}}\right)^{\frac{2}{\alpha }},\gamma =\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{(2n+1)^{1+\frac{2}{\alpha }}}.$$
(2)
Obviously, if the 1d system is isotropic, the MFPT for the carrier travelling from the right to the left end $`t_{RL}`$ is the same as $`t_{LR}`$.

If the heat bath is a stochastic kernel of Gaussian type, namely, the probability distribution of velocities is $`p(v,T)=4\pi v^2\mathrm{exp}(v^2/2T)/(2\pi T)^{3/2}`$, the MFPT becomes,

$$t_{LR}=\frac{16T\gamma }{\alpha (2\pi T)^{\frac{3}{2}}}\left(\frac{2L}{\pi \sqrt{D}}\right)^{\frac{2}{\alpha }}.$$
(3)
We define the heat current as the energy exchange between two heat baths in unit time. Thus the current induced by a carrier ($`m=1`$) with velocity $`v`$ moves from left to right and comes back is:

$$j=\frac{_0^{\mathrm{}}\frac{v^2}{2}\left(p(v,T_L)p(v,T_R)\right)}{t_{LR}+t_{RL}}=\frac{T_LT_R}{2t_{LR}},$$
(4)
If the temperature difference between the two baths is sufficient small so that $`T=(T_RT_L)/L`$, then the thermal conductivity, $`\kappa =Lj/T`$, is

$$\kappa =cL^\beta ,\beta =22/\alpha ,$$
(5)
and the constant, $`c=3\pi \sqrt{2\pi }\alpha (\pi \sqrt{D}/2)^{2/\alpha }\sqrt{T}/(32\gamma )`$.

Eq. (5) is the central result of the paper. It connects heat conduction and diffusion quantitatively. The main conclusion is that an anomalous diffusion indicates an anomalous heat conduction with a divergent (convergent) thermal conductivity. More precisely, our result tells us that: a ballistic motion means thermal conductivity proportional to the system size $`L`$, a normal diffusion means a normal heat conduction obeying the Fourier law, a superdiffusion means a divergent thermal conductivity, a subdiffusion means a zero thermal conductivity in the thermodynamic limit. In the following, we compare our results with the existing analytical and numerical results.

A ballistic motion, $`\alpha =2`$, leads to a divergent thermal conductivity $`\kappa L`$. The only existing analytical result is heat conduction in a 1D harmonic lattice. It is known that heat is transported by phonons in lattice model. Because there are no resistance and umklapp process, the phonons transport ballistically in harmonic lattice model, thus $`\alpha =2`$. From our formula (5), the thermal conductivity in the 1D harmonic lattice diverges as $`L^\beta `$ with $`\beta =1`$, this is exactly what was shown by Lebowitz et al (“$``$” in Fig. 1).

A normal diffusion, $`\alpha =1`$, means that the thermal conductivity is a size independent constant, $`\beta =0`$, i.e. the heat conduction obeys the Fourier law. For example in the 1D Frenkel-Kontorova model, in a certain range of parameter such as $`T/K<<1`$, the phonons transport diffusively, thus the thermal conductivity is finite and independent of the system size $`L`$. The disordered FPU model also has a finite thermal conductivity due to the random walk like scattering process in the chain. Other 1D models showing normal diffusion and normal thermal conduction are: the 1D Lorentz gas channel, the 1D disordered Ehrenfest gas channel, the 1D irrational triangle channel, the alternative mass hard-core potential model, and some 1D polygonal billiard channels with certain rational triangles. “$``$” in Fig. 1 represents all models with normal diffusion.

A superdiffusion, $`1<\alpha <2`$, implies an anomalous heat conduction with a divergent thermal conductivity $`L^\beta `$. The exponent $`0<\beta <1`$ differs from model to model. Here we take the billiard models as our examples because they are very clean, and both the diffusion and thermal conductivity in these models can be calculated very accurately. The first example is the 1D Ehrenfest gas channel in which the scattering obstacles are isosceles right triangles periodically post along the channel. In this model one has $`\alpha =1.672`$. From our analytical result (5), the thermal conductivity should diverges as $`L^\beta `$ with $`\beta =22/\alpha =0.804`$ which agrees with the result from simulation of heat conduction $`\beta =0.814`$ (“$``$” in Fig. 1). The second example is the 1D channel with triangles whose inter angles are rational multiples of $`\pi `$. This model shows a superdifussion with $`\alpha =1.178`$. The divergent exponent of thermal conductivity is $`\beta =0.302`$. This exponent is slightly larger than the one obtained from thermal conductivity simulation $`\beta =0.22`$ (“$``$” in Fig. 1). This deviation is due to the finite size effect in the heat conduction simulation.

A subdiffusion, $`\alpha <1`$, results in an anomalous heat conduction with a convergent thermal conductivity, i.e. $`\kappa L^\beta `$, with $`\beta <0`$. This is an interesting result implying that the system becomes a thermal insulator in the thermodynamic limit $`L\mathrm{}`$. Although there are many examples showing subdiffusion , a systematic study on the heat conduction in such kind of systems is still lacking. The only existing example is the heat transport in a polygonal billiard which supports our conclusion. Most recently, Alonso et al show that in a very special configuration, $`\alpha =0.86`$, and the thermal conductivity goes as $`\kappa L^{0.63}`$ (“$`\mathrm{}`$” in Fig. 1). As $`L`$ goes to infinity, the thermal conductivity goes to zero. According to our formula (5), if $`\alpha =0.86`$, $`\beta =0.33`$ which is larger than the one obtained by Alonso et al. This is not a surprise, because the channel length in their study of thermal conductivity is too small ($`L40`$). If the channel is longer, the value of $`\beta `$ will become much more closer to our theoretical estimation ($`\beta =0.33`$). To demonstrate this, we extend the thermal conductivity simulation from $`L[1,40]`$ used by Alonso et al  to $`L[40,80]`$, and find that $`\beta =0.48`$ (“$``$” in Fig. 1) which is more closer to $`\beta =0.33`$ than the one obtained by Alonso et al. If $`L\mathrm{}`$, one can expect $`\beta `$ goes to $`0.33`$.

All numerical results are summarized and represented in Fig 1, where we draw $`\beta `$ versus $`\alpha `$, and compared with Eq. (5). As is shown that, Eq. (5) is exact for both normal diffusion and the ballistic motion. The agreement with most existing numerical data is good. However, discrepancies remain for some models mainly due to the limited numerical simulations. The best data close to curve $`\beta =22/\alpha `$ is the simulation from 1D Ehrenfest gas channel. This is because the channel length used in the simulation is the longest one ($`L10^3`$) among all the available data.

In summary, we have established a connection between anomalous heat conduction and anomalous diffusion in 1D systems. Our central result Eq. (5) includes all possible cases observed in different classes of 1D models, ranging from subdiffsion, normal diffusion, and superdiffusion to ballistic transport. Several conclusions can be drawn: (1) A normal diffusion leads to a normal heat conduction obeying the Fourier law. (2) A ballistic transport leads to an anomalous heat conduction with a divergent thermal conductivity $`\kappa L`$. (3) A superdiffusion leads to an anomalous heat conduction with a divergent thermal conductivity in thermodynamic limit. (4) More importantly, our result predicts that a subdiffusion system will be a thermal insulator. Existing numerical data support our results.

We should mention that the subdiffusion process has been observed in many real physical systems such as highly ramified media in porous systems, percolation clusters, exact fractals, the motion of a bead in a polymer network, charge carrier transport in amorphous semiconductors. Any numerical simulation and/or real experimental measurement of thermal conductivity in these systems will be very interesting and will allow one to test the theory given in this Letter. More importantly, it will have a wide application in designing novel thermal devices.

We would like to thank D Alonso for providing us reference before publication, and G Casati for using discussions. This work was supported in part by Academic Research Fund of NUS. JW is supported by DSTA Singapore under Project Agreement POD0001821. This paper has been reported at the international workshop and seminar ”Microscopic Chaos and Transport in Many-Particle Systems”, at the MPIPKS Dresden from August 5 to August 25, 2002.