text:
# Comment on “Measuring the Orbital Angular Momentum of a Single Photon”

## Abstract

Optical modes with different orbital angular momentums (OAMs) per photon may be sorted by Mach-Zehnder interferometers incorporated with beam rotators, without resorting to OAM mode converters.

Presented in a recent letter Leach02  is an ingenious method of sorting spatial modes of photons with different orbital angular momentum (OAM), which closely resembles, and in a sense complements a scheme of analyzing optical beams with rectangular symmetry Xue01 . The method in Leach02  employs special Mach-Zehnder (MZ) interferometers as binary branching devices, each of which divides a set of input modes with $`lk(\mathrm{mod}2^n)`$ into two groups with $`lk(\mathrm{mod}2^{n+1})`$ and $`l2^n+k(\mathrm{mod}2^{n+1})`$ respectively, where the integer $`l`$ denotes the single-photon OAM in units of $`\mathrm{}`$, $`k`$ and $`n`$ are fixed, non-negative integers specific to the individual interferometer. In the Letter, the same condition is assumed implicitly for all the MZ interferometers, that one arm maintains zero phase (or integral multiples of $`2\pi `$) for all the modes, while the other arm rotates the beam by an angle $`\alpha `$ so to induce an OAM-dependent phase shift $`l\alpha `$. The assumption proves rather restrictive and responsible for the necessity of OAM mode converters, which result in optical loss and increase significantly the complexity of the optical setup. The loss of light could impose a serious limitation to quantum optical experiments involving single photons.

The use of OAM mode converters may be avoided by relaxing the aforementioned restriction on the MZ interferometers, namely, by incorporating an adjustable phase shifter in a modified MZ interferometer as depicted in Fig.1. In practice, there is always a mechanism of phase adjustment in setting up an optical interferometer. An example implementation is a thin glass film making an adjustable angle to the beam axis. When OAM modes with $`lk(\mathrm{mod}2^n)`$ are input into the modified MZ interferometer, the beam is rotated by $`\alpha =\pi /2^n`$ in the upper arm such that a mode with OAM $`l`$ acquires a phase $`l\alpha =l\pi /2^n`$, while the phase shifter in the lower arm is tuned to induce a fixed phase shift $`k\alpha =k\pi /2^n`$ to all the OAM modes. It is easily seen that the modified MZ interferometer segregates the OAM modes into two groups with $`lk(\mathrm{mod}2^{n+1})`$ and $`l2^n+k(\mathrm{mod}2^{n+1})`$ respectively. Such modified MZ interferometers may be used as binary branching devices to construct an OAM mode sorter, without resorting to OAM mode converters.

The OAM mode sorter is in striking similarity to the Hermite-Gaussian (HG) mode analyzer using MZ interferometers incorporated with fractional Fourier transformers (FRFTs) Xue01 . FRFT-incorporated MZ interferometers can even be used after an OAM mode sorter to lift the mode degeneracy due to the radial degree of freedom, so that orthogonal Laguerre-Gaussian (LG) modes can be sorted completely. Moreover, a complete HG mode sorter may be utilized just as a complete LG mode analyzer with the help of HG $``$ LG mode converters Beijersbergen93 . Finally, it may be noted that the increased information capacity using higher order OAM Leach02 ; Molina-Terriza02  or HG modes is just the result of an enlarged spatial channel Miller98 . To achieve the higher capacity, the spatial channel needs to support optical beams with larger sizes than the fundamental mode.


format: text

text:
# X-ray Spectroscopy of the Accreting Millisecond Pulsar XTE J0929–314 in Outburst

## 1 Introduction

Millisecond pulsars (MSPs) have long been considered one of the possible endpoints of low-mass X-ray binary (LMXB) evolution. The neutron star (NS) is thought to be spun-up to a millisecond period by accretion from its low-mass companion. After the accretion phase has ended, the NS may turn on as a radio MSP. In the last five years, this theory has been confirmed with the identification of three accretion-powered MSPs, SAX J1808.4$``$3658, XTE J1751$``$305 and XTE J0929$``$314 (Wijnands & van der Klis, 1998; Chakrabarty & Morgan, 1998; Markwardt et al., 2002; Galloway et al., 2002).

Interestingly, all three accretion-powered MSPs are in short period binaries, with the two most recently discovered having binary periods $``$43 min (Markwardt et al., 2002; Galloway et al., 2002). These short periods place the MSPs XTE J1751$``$305 and XTE J0929$``$314 in the class of ultracompact binaries, defined as having orbital periods $``$80 min. Ultracompact binaries require hydrogen-deficient or degenerate donors (Nelson, Rappaport, & Joss, 1986). Recently, O and Ne emission and absorption features were discovered in two known and three suspected ultracompact systems (Schulz et al., 2001; Juett, Psaltis, & Chakrabarty, 2001; Juett & Chakrabarty, 2002). These results led the authors to conclude that the donor stars are degenerate C-O WDs. On the other hand, the XMM-Newton spectrum of XTE J1751$``$305 did not show emission or absorption features (Miller et al., 2002). In addition, no unusual abundances were required to fit the neutral absorption edges. This result is consistent with the suggestion of Markwardt et al. (2002) that the donor in XTE J1751$``$305 is a He WD (see also, Bildsten, 2002).

XTE J0929$``$314 was discovered in April 2002 by the All Sky Monitor onboard the Rossi X-ray Timing Explorer (RXTE; Remillard, 2002). Further RXTE observations detected 185 Hz pulsations modulated by a 43.6-min binary orbit (Remillard, Swank, & Strohmayer, 2002; Galloway et al., 2002). Radio and optical counterparts with positions consistent with the X-ray position were also detected (Rupen, Dhawan, & Mioduszewski, 2002; Greenhill, Giles, & Hill, 2002; Cacella, 2002). An optical spectrum of XTE J0929$``$314 revealed emission lines from C III/N III $`\lambda `$4640-4650 and H$`\alpha `$ $`\lambda `$6563 (Castro-Tirado et al., 2002). Given the optical detection of emission lines, and the ultracompact nature and high-Galactic latitude of the source, XTE J0929$``$314 is an ideal candidate for a high-resolution X-ray spectroscopic study to search for emission and absorption features similar to those seen in other ultracompact systems. In this letter, we present results from a Director’s Discretionary Time observation of XTE J0929$``$314 with the Chandra X-ray Observatory, as well as spectral results from the RXTE pointed observations throughout the outburst.

## 2 Observation and Data Reduction

The 2002 outburst of XTE J0929$``$314 began around MJD 52370 and peaked at $`31`$ mCrab (2–10 keV) on MJD 52394 (Galloway et al., 2002). The pointed RXTE observations commenced 2 d later, and proceeded once every few days through MJD 52456. The source dropped below the $`3\sigma `$ detection limit of $`7.5\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ (2–10 keV; equivalent to 0.32 mCrab) after MJD 52443.

XTE J0929$``$314 was observed throughout its 2002 outburst by the two pointed instruments onboard RXTE. The Proportional Counter Array (PCA; Jahoda et al., 1996) consists of five gas-filled proportional counter units (PCUs) with a total effective area of $`6000\mathrm{cm}^2`$, sensitive to X-ray photons in the 2.5–60 keV range. The two instrument clusters comprising the High-Energy X-ray Timing Experiment (HEXTE; Gruber et al., 1996) present an effective area of $`1600\mathrm{cm}^2`$ to photons in the energy range 15–250 keV. We extracted average spectra from standard instrument mode data (“Standard-2” and “Archive” for PCA and HEXTE, respectively) from each PCU/cluster for each observation. Since the instrumental gain is known to vary between PCUs, and also with time, we generated a separate response matrix for each PCU and each observation using pcarsp version 8.0, supplied with LHEASOFT version 5.2. We estimated background spectra using the “combined” gain epoch 5 (beginning 2000 May 13) faint-source model as input to pcabackest version 3.0. Over the course of the outburst, we collected 38 observations of XTE J0929$``$314 totaling 123.5 ks.

We also observed XTE J0929$``$314 with Chandra on 2002 May 15 (MJD 52409) for 18 ks using the Low Energy Transmission Grating Spectrometer (LETGS) and the Advanced CCD Imaging Spectrometer (ACIS; Weisskopf et al., 2002). The LETGS spectra are imaged by ACIS, an array of six CCD detectors. The LETGS/ACIS combination provides both an undispersed (zeroth order) image and dispersed spectra from the grating with a first order wavelength range of 1.4–63 Å (0.2–8.9 keV) and a spectral resolution of $`\mathrm{\Delta }\lambda =`$ 0.05 Å. The various orders overlap and are sorted using the intrinsic energy resolution of the ACIS CCDs. The observation used a Y-offset of $`+`$1$`\stackrel{}{\mathrm{.}}`$5 in order to place the O-K absorption edge on the back-side illuminated S3 CCD, which has suffered less degradation than the front-side illuminated CCDs.

Using tgdetect, we determined the zeroth order source position of XTE J0929$``$314 to be: R.A.=$`09^\mathrm{h}29^\mathrm{m}20\stackrel{\mathrm{s}}{\mathrm{.}}15`$ and Dec=$`31^{}23\mathrm{}04\stackrel{}{\mathrm{.}}3`$, equinox J2000.0 (90% confidence error of 0$`\stackrel{}{\mathrm{.}}`$6111See http://asc.harvard.edu/cal/ASPECT/celmon/index.html). The Chandra position is consistent with both the optical and radio counterpart positions (Rupen et al., 2002; Greenhill et al., 2002; Cacella, 2002). The first order dispersed spectrum of XTE J0929$``$314 had an average count rate of 7.64$`\pm `$0.02 counts s-1. We examined the total count rate, as well as the count rates in two different energy ranges, to check for changes in the spectral state. We found no evidence for any change of spectral state during the Chandra observation.

The “level 1” event file was processed using the CIAO v2.2 data analysis package222http://asc.harvard.edu/ciao/. The standard CIAO spectral reduction procedure was performed. We filtered the event file retaining those events tagged as afterglow events by the acis\_detect\_afterglow tool. Since order-sorting of grating spectra provides efficient rejection of background events, the afterglow detection tool is not necessary to detect cosmic ray afterglow events. No features were found that might be attributable to afterglow events.

For bright sources, pileup can be a problem for CCD detectors (see, e.g., Davis, 2002). The zeroth order Chandra spectra of XTE J0929$``$314 was heavily affected by pileup and was not used in this analysis. In addition, it was found that the first order spectrum of XTE J0929$``$314 suffered from pileup in the range 5–12 Å (1–2.5 keV). In order to use the grating pileup kernel in ISIS, we created detector response files (ARFs) for each chip using the CIAO tool mkgarf and then combined them to create $`+1`$ and $`1`$ order ARFs using a custom tool developed by J. Davis of the HETGS instrument team. (This tool is similar to the standard CIAO tool but correctly calculates the fractional exposure at each response bin. This functionality has been added to the standard tools in CIAO v2.3.) The grating pileup kernel models the effect of pileup on grating spectra, similar to the pileup model available for CCD spectra in ISIS and XSPEC (see, Davis, 2002, for a discussion of grating pileup modeling). In addition, the data and responses were rebinned to 0.083 Å to reflect the size of an ACIS event detection cell (3 CCD pixels). Background spectra were extracted from the standard LETG background regions. Spectral analysis of the Chandra observation of XTE J0929$``$314 was performed using ISIS (Houck & DeNicola, 2000).

## 3 Search for Orbital Modulation

The presence (or absence) of orbital modulation of the X-ray flux can help constrain the inclination of the binary, and in turn the mass of the companion. Galloway et al. (2002) determined the orbital period of XTE J0929$``$314 to be 43.6 min from the Doppler modulation of the pulse arrival times. We searched both the RXTE and Chandra data for orbital modulation of the X-ray flux.

For the RXTE analysis, we calculated background-subtracted, solar-system barycenter corrected, 16-s binned lightcurves for each observation from photons in the energy range 3–13 keV. This range was selected to maximize the signal-to-noise ratio, and is identical to that used for the pulse timing analysis of Galloway et al. (2002). For each observation, we created a phase folded lightcurve of 8 bins using the measured orbital parameters, with phase zero set to the epoch of $`90^{}`$ mean longitude, $`T_{\pi /2}`$. The length of each observation was short compared to the decay timescale of the X-ray flux of XTE J0929$``$314, so it was not necessary to correct for the decline in X-ray flux during a single observation. Each lightcurve was then fit with a sine curve with variable amplitude and phase. From these results, we place a 3$`\sigma `$ upper limit on the fractional rms modulation of 1.1%.

We also searched the Chandra data for orbital modulation of the X-ray flux. For the Chandra analysis, we first barycentered and randomized the first order grating events. Randomizing of the event arrival times consists of adding a random quantity uniformly distributed between 0–3.2 s, in order to avoid aliasing caused by the readout time. We calculated the phase of each first order event using the orbital parameters from Galloway et al. (2002) and created an 8 bin phase-folded lightcurve. We found an upper limit of 2.3% rms, consistent with the RXTE results.

Since the Chandra data is not subject to time gaps due to spacecraft motion, we also calculated a Fourier transform of the first order events which had been made into a 10-s binned lightcurve. We searched for modulations of the X-ray flux with frequencies between 6$`\times 10^5`$ and 5$`\times 10^2`$ Hz. We found no evidence for any periodic modulation, with a 90%-confidence upper limit of 1.4% for the fractional rms amplitude over the frequency range. The 90%-confidence upper limit for the fractional rms amplitude of a signal at the orbital period is 0.4%, well within the upper limits given by the phase folding analyses.

## 4 Spectral Analysis

### 4.1 RXTE Spectral Analysis

We fitted the combined 2–25 keV PCA and 15–50 keV HEXTE spectra with a number of models which typically give good fits for other X-ray pulsars, including broken and cutoff power laws and analytical Comptonization approximations. There was no evidence for a high-energy spectral cutoff within the energy range in which the source was detected (typically $`50`$ keV). We did however measure significant residuals below 10 keV. These residuals could be minimized by adopting a broken power law model where the spectral index decreased by $`10`$% (i.e. the spectrum became harder) above 7 keV, or by adding a blackbody component with $`kT_{\mathrm{bb}}0.5`$ keV to a power law. While both these models gave fits of similar quality, we rejected the former as unphysical and instead adopted the latter for the remainder of the RXTE fits.

We also included a multiplicative component to take into account the attenuation by intervening (neutral) material with cosmic abundances. The equivalent hydrogen column density $`N_\mathrm{H}`$ exhibited no significant variations between observations, and in the mean was $`10^{21}\mathrm{cm}^2`$, consistent with the line-of-sight values interpolated from dust and H I maps ($`7.6\times 10^{20}`$ and $`10^{21}\mathrm{cm}^2`$ respectively; Schlegel, Finkbeiner, & Davis, 1998; Dickey & Lockman, 1990). For consistency with the fits to Chandra data, we froze the $`N_\mathrm{H}`$ value in subsequent fits at $`7.6\times 10^{20}\mathrm{cm}^2`$. The mean value of the reduced-$`\chi ^2`$ ($`=\chi _\nu ^2`$) for power-law $`+`$ blackbody fits to spectra from all the observations was 1.04. While the worst-fitting observation gave a $`\chi _\nu ^2=2.44`$, the largest residuals arose from variations between PCUs (particularly at the lower energy bound) rather than any broad trend with energy. By excluding these lowest channels from the fit, and assuming a systematic error of 1%, we were able to reduce the fit statistic for that observation to $`\chi _\nu ^2=1.30`$.

The weighted mean of the blackbody temperature, $`kT_{\mathrm{bb}}`$, was $`0.66\pm 0.06`$ keV, while the normalization $`\sqrt{K_{\mathrm{bb}}}=R_{\mathrm{bb},\mathrm{km}}/d_{10\mathrm{k}\mathrm{p}\mathrm{c}}=5.9\pm 1.3`$ (all errors are quoted at 90%-confidence unless otherwise noted). Both parameters showed variations during the outburst; however, they were significantly anticorrelated, so that the estimated bolometric blackbody flux ($`K_{\mathrm{bb}}T_{\mathrm{bb}}^4`$) remained constant to within our measurement uncertainties ($`\chi _\nu ^2=1.06`$, see Figure 1). In the first observation, the blackbody component comprises only 3% of the flux in the 2–60 keV band, whereas in the last observation where significant blackbody and power-law components were detected, the blackbody accounted for 16% of the flux.

We calculated a weighted mean photon index of $`\mathrm{\Gamma }=1.76\pm 0.03`$. The principal source of flux variation was the normalization of the power-law component (see Figure 1). While all parameters exhibited variations throughout the outburst, the power-law normalization was the only parameter to show a significant correlation with the integrated 2–10 keV flux (Spearman’s rank correlation coefficient $`\rho =0.908`$, equivalent to 5.4$`\sigma `$; Press et al., 1992). A power-law component was not significantly detected in the last observation on MJD 52443. A blackbody-only fit for that day exhibits a significantly (6.6$`\sigma `$) lower integrated blackbody flux than the weighted mean level for the previous observations. We did not detect XTE J0929$``$314 in subsequent observations (3$`\sigma `$ upper limit on the 2–10 keV flux of $`7.5\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$).

### 4.2 Chandra Spectral Analysis

We fit the Chandra spectrum with a power-law $`+`$ blackbody model including absorption. The $`+`$1$`\stackrel{}{\mathrm{.}}`$5 Y-offset in the observational setup allows for the entire 1.5–25.3 Å (0.5–8.3 keV) $`+1`$ order spectrum to be imaged on the S3 CCD. This CCD is backside-illuminated and has suffered less degradation than the front-side illuminated CCDs. Because of this and the lack of counts above 25 Å, we performed our spectral fits on only the +1 order spectrum in the range 1.5–25.3 Å.

Of the three known X-ray MSPs, XTE J0929$``$314 has the highest Galactic latitude and, probably the smallest interstellar column density along the line of sight to the source. Our Chandra observation allows for a direct determination of the strength of the absorption edges, in particular the O-K edge and the Fe-L triplet. To account for absorption and measure the optical depths of the edges, we used the tbvarabs absorption model (see, Wilms, Allen, & McCray, 2000) with O and Fe abundances set to zero and the $`N_\mathrm{H}`$ fixed to 7.6$`\times 10^{20}`$ cm-2, the hydrogen column from dust maps (Schlegel et al., 1998). The O edge was fit with an edge model, while the Fe-L edges were fit using a custom model that employs the optical constant measurements of Kortright & Kim (2000). In addition, the C abundance in tbvarabs was allowed to vary in order to account for the known instrumental contamination of the ACIS CCDs333For more information see, http://cxc.harvard.edu/cal/Links/Acis/acis/Cal\_prods/qeDeg/index.html (Plucinsky et al., 2002). We also included Gaussian lines to fit the interstellar atomic O absorption line at 23.5 Å and the absorption feature

at 23.36 Å, attributed to a 1s-2p transition in Fe2O3 (see, Schulz et al., 2002).

We found a best-fit photon index of 1.55$`\pm `$0.03, with a normalization of (4.04$`\pm `$0.15)$`\times 10^2`$ photons keV-1 cm-2 s-1 at 1 keV. The Chandra best-fit blackbody temperature was $`kT_{\mathrm{bb}}=0.65\pm 0.03`$ keV, with $`R_{\mathrm{bb},\mathrm{km}}/d_{10\mathrm{k}\mathrm{p}\mathrm{c}}=7.6\pm 0.8`$. Our best-fit model had a $`\chi _\nu ^2=1.001`$ and an unabsorbed 2–10 keV flux of 2.7$`\times 10^{10}`$ erg cm-2 s-1.

The best-fit optical depth for O was 0.34$`\pm `$0.07, which translates to an equivalent hydrogen column density of (1.2$`\pm `$0.3)$`\times 10^{21}`$ cm-2, using the cross-section of Henke, Gullikson, & Davis (1993) and the O ISM abundance of Wilms et al. (2000). Although this is somewhat higher than the measured hydrogen column from dust maps, an instrumental contribution to the O edge has been reported, with an optical depth of $``$0.10 (Plucinsky et al., 2002). This contribution would add $`N_\mathrm{H}0.36\times 10^{21}`$ cm-2 to the expected value, making the Chandra O edge measurement compatible with the dust map measurement. The best-fit line wavelength and equivalent width (EW) of the atomic O absorption line were 23.52$`\pm `$0.05 Å and $`0.053_{0.03}^{+0.005}`$ Å. The Fe2O3 absorption line and Fe-L edges were not significantly detected, with an upper limit on the Fe column density of $`3\times 10^{16}`$ cm-2, equivalent to an $`N_\mathrm{H}<1.1\times 10^{21}`$ cm-2. The best-fit C abundance was $`7.8_{0.7}^{+1.2}`$ times the interstellar C/H ratio of Wilms et al. (2000). This is slightly lower than the expected instrumental $`+`$ interstellar contribution of 9.9$`\pm `$0.5 times the interstellar ratio. The instrumental edge depths for C and O were estimated using the recent calibration results presented at the Chandra X-ray Center website55footnotemark: 5 (see also, Plucinsky et al., 2002).

Besides the interstellar O absorption line, there were no prominent emission or absorption lines in the high-resolution spectrum of XTE J0929$``$314 (see Figure 2). We performed a careful search of the Chandra spectral residuals to place limits on the presence of any spectral features. Gaussian models with fixed FWHM $`=2000`$ km s-1 were fit at each point in the wavelength range 1.5–25.3 Å. From this, we can place a 3$`\sigma `$ upper limit on the EW of any line feature, either emission or absorption. The EW limit increases approximately linearly with wavelength, varying from 0.007 Å at 1.5 Å to 0.12 Å at 24 Å.

Our measured C abundance depends on the assumed continuum model since the edge depth can not be directly measured (the C edge at 43 Å is outside our bandpass). In order to test if the C abundance was the cause of the difference in the photon index between the Chandra and RXTE fits, we refit the Chandra spectrum with a fixed C abundance of 9.9 times the interstellar C/H ratio. With the C abundance fixed, the best-fit photon index was 1.62$`\pm `$0.03. This is still significantly different from the best-fit photon index found in the RXTE spectra. The best-fit power-law normalization increased to (4.50$`\pm `$0.07)$`\times 10^2`$ photons keV-1 cm-2 s-1. The other fit parameters were consistent with the previous fit within errors. We also considered the possibility that the pileup model was giving rise to the lower photon index in the Chandra spectral fits. When the Chandra data is fit without the pileup model, and excluding the piled-up region (1–3 keV) of the spectrum, we find spectral parameters consistent with the results found using the pileup model.

## 5 Combined Chandra and RXTE Spectral Fits

To investigate the difference in the spectral results between Chandra and RXTE, we performed combined spectral fits with the Chandra and RXTE data. A single RXTE observation (70096-03-05-00), taken within a few hours of the Chandra data, was used to represent the RXTE spectrum of XTE J0929$``$314. To simplify fitting, the piled-up region (1-3 keV) of the Chandra data was excluded and the pileup kernel was not used. The combined data were fit with a blackbody $`+`$ either a power-law, broken power-law or the comptt model which describes Comptonization of soft photons by a hot plasma (Titarchuk, 1994). The broken power-law and Comptonization models allow for turnover at low energies which could explain the difference in the Chandra and RXTE spectral results. The absorption was described by the tbvarabs model with the C abundance allowed to vary for the Chandra observation only, while best-fit edge depths from the Chandra analysis were used for the O and Fe edges. We also included a Gaussian line, fixed to the best-fit Chandra value, to model the atomic O absorption at 23.5 Å. Multiplicative constants were included to compensate for any instrumental normalization differences. We find that the normalizations of the PCA spectra are $``$20% larger than for Chandra, while the HEXTE spectra are $``$20% lower. Similar normalization differences have been noted before (see, e.g., Kuulkers et al., 2002, and references therein).

As a check on the instrumental calibrations, we directly compared the Chandra and RXTE data in the same energy range by fitting the high energy (3–8 keV) Chandra data, in conjunction with the RXTE data, to an absorbed power-law $`+`$ blackbody model with absorption fixed to $`N_\mathrm{H}=7.6\times 10^{20}`$ cm-2. The fits are consistent with the results of the RXTE data, $`\mathrm{\Gamma }=1.74\pm 0.02`$, $`kT_{\mathrm{bb}}=0.66\pm 0.03`$ keV, and $`R_{\mathrm{bb},\mathrm{km}}/d_{10\mathrm{k}\mathrm{p}\mathrm{c}}=5.9\pm 0.7`$, with $`\chi _\nu ^2=1.03`$.

We then fit the full Chandra and RXTE spectra with a power-law $`+`$ blackbody model. When the edge parameters were fixed to the Chandra best-fit values, the photon index had a best-fit value of 1.628$`\pm `$0.009 and with a $`\chi _\nu ^2=1.20`$ (see Table 1). We note that the $`\chi _\nu ^2`$ of the power-law $`+`$ blackbody fit to the RXTE data alone was 1.06. If we fix the C abundance to the expected value of 9.9 times the interstellar C/H ratio, the photon index increases to 1.656$`\pm `$0.009 ($`\chi _\nu ^2=1.14`$). If we instead allow the C abundance to vary, the best-fit photon index increases to 1.72$`\pm `$0.02 ($`\chi _\nu ^2=1.09`$), giving a photon index consistent with the RXTE spectral results, but with a C abundance of 14.1$`\pm 1.2`$ times the interstellar C/H ratio. This value is inconsistent with the estimate of the C abundance from instrument calibration work (see §4.2). As can be seen, the continuum spectral parameters are highly dependent on the assumed absorption.

We performed the same fits using the comptt in place of the power-law component. The Comptonization model provides for a low-energy turnover in the spectrum. The electron temperature was fixed to 50 keV, the upper limit of the energy range. In addition, we used a spherical geometry in order to compare with the spectral fit of SAX J1808.4$``$3658 (Titarchuk, Cui, Wood, 2002). The comptt model provides a better fit to the data compared with the power-law ($`\chi _\nu ^2=1.09`$–1.07, see Table 1). We performed an $`F`$-test to formally compare the two models and found that in all cases the Comptonization model was better with a significance $`>99`$%. For the two fixed values of the C abundance, the results of the comptt $`+`$ blackbody fits are consistent within error with a seed photon temperature, $`kT_00.13`$ keV and optical depth, $`\tau _p2.7`$. When the C abundance is allowed to vary, the fit improves but becomes inconsistent with the other results, which is probably due to degeneracy between the parameters.

In the best-fit comptt $`+`$ blackbody model, the blackbody component has parameter values of $`kT_{\mathrm{bb}}=0.71`$ keV and $`R_{\mathrm{bb},\mathrm{km}}/d_{10\mathrm{k}\mathrm{p}\mathrm{c}}=5.5`$. These values are consistent with the power-law fits both to the combined data and with the blackbody parameters from each instrument fit, indicating that the blackbody temperature and flux are independent of the exact continuum model used. We find $`kT_0<kT_{\mathrm{bb}}`$, similar to the fit of SAX J1808.4$``$3658 using the same model (Titarchuk et al., 2002). We note that the comptt model is in itself an approximation to a full treatment of Comptonization, in particular at soft energies. If we assume that the blackbody provides the input photons to the Comptonization component, then we would expect $`kT_0=kT_{\mathrm{bb}}`$. We fit the data with the constraint $`kT_0=kT_{\mathrm{bb}}`$, but the fit was unreasonable with $`\chi _\nu ^2>2.0`$. Interestingly, Gierliński, Done, & Barret (2002) used a different Comptonization model and found $`kT_0=kT_{\mathrm{bb}}`$ in spectral fits of SAX J1808.4$``$3658. We employed the same model to fit XTE J0929$``$314 (compPS; Poutanen & Svensson, 1996), but found similar values for the input photon temperature and optical depth as found using the comptt model.

We also fit the data with a broken power-law $`+`$ blackbody model to test for low-energy turnover without relying on the assumptions inherent in the Comptonization model. We find that the data is also well fit ($`\chi _\nu ^2=1.08`$) by a broken power-law $`+`$ blackbody model with $`\mathrm{\Gamma }_1=1.26`$–1.56, $`\mathrm{\Gamma }_2=1.75`$, and $`E_{\mathrm{break}}=1.4`$–4.4 keV, and blackbody parameters consistent with previous fits (see Table 1).

## 6 Discussion

We have found that the spectrum of XTE J0929$``$314 is well described by a power-law $`+`$ blackbody model over limited energy ranges, with interstellar absorption consistent with $`N_\mathrm{H}=7.6\times 10^{20}`$ cm-2, the expected hydrogen column density along the line of sight to the source. The Chandra spectrum shows no prominent emission or absorption lines in the spectrum (similar to XTE J1751$``$305; Miller et al., 2002), with an EW limit which varies from $`<0.007`$ Å at 1.5 Å to $`<0.12`$ Å at 24 Å. No orbital modulation of the X-ray flux was detected in either the Chandra or RXTE data. From the lack of X-ray eclipses, we set an upper limit on the binary inclination of $`i85^{}`$ for a Roche-lobe–filling companion. In addition, the absence of dips in the X-ray lightcurve suggests that the source is most likely at an inclination of $`i60^{}70^{}`$ (see, e.g., Frank, King, & Lasota, 1987; White, Nagase, & Parmar, 1995).

The blackbody spectral component in LMXBs is generally attributed to the neutron star, possibly from a “hot spot” where the accretion column meets the neutron star surface. From our fits, we find a blackbody radius of $`R_{\mathrm{bb},\mathrm{km}}=(5.37.7)d_{10\mathrm{k}\mathrm{p}\mathrm{c}}`$. Given the lower limit on the distance of 5 kpc (Galloway et al., 2002), we find a lower limit on the blackbody radius of $`R_{\mathrm{bb}}>2.7`$ km. We note that this lower limit does not include a correction for the conversion between the color temperature and the effective temperature of the neutron star atmosphere (see, Lewin, van Paradijs, & Taam, 1993, for a discussion). This correction could increase the inferred radius by a factor of $`2`$. But even after this correction, the inferred radius is not consistent with canonical NS models for a distance of 5 kpc suggesting a “hot spot” emission region.

During the outburst, we found that the blackbody flux remained constant, while the power-law flux declined. In addition, when the combined Chandra and RXTE spectrum is fit with a Comptonization model, we find that the blackbody temperature, $`kT_{\mathrm{bb}}`$, is significantly greater than the seed photon temperature of the Comptonization model, $`kT_0`$. This evidence suggests that the power-law and blackbody components are associated with different emission processes.

We compared our results with two physical models that were developed to explain the spectral properties of the accretion-powered MSP SAX J1808.4$``$3658 (Gierliński et al., 2002; Titarchuk et al., 2002). Both models were formulated using data from the 1998 outburst of SAX J1808.4$``$3658, but come to very different conclusions regarding the relationship between the power-law or Comptonization component and the blackbody. Gierliński et al. (2002) assumed that the blackbody was the source of the seed photons (i.e, $`kT_0=kT_{\mathrm{bb}}`$) and proposed that the Comptonization component comes from the shock heated accretion column with soft photon input from the “hot spot” blackbody component. Such a model implies that the blackbody and power-law/Comptonization component fluxes should be correlated, which was the case for SAX J1808.4$``$3658 but was not seen for XTE J0929$``$314.

On the other hand, Titarchuk et al. (2002) found that $`kT_0<<kT_{\mathrm{bb}}`$ and suggested that the accretion disk provided the majority of the soft photon input to a spherical Compton cloud, of radius $`R60`$ km, from which the Comptonization component originates. The Titarchuk et al. (2002) model allows for the the blackbody and power-law/Comptonization component to vary independently. Titarchuk et al. (2002) also suggested that the low ($`\tau _04`$) optical depth of the Comptonizing region found for SAX J1808.4$``$3658 allows for the pulsations of the NS to be detected, while the higher optical depths ($`\tau _0>5`$) reported in more luminous LMXBs suppress the pulsation amplitudes. We find a best-fit $`\tau _02.7`$ for XTE J0929$``$314, consistent with this explanation. However, we note that Miller et al. (2002) did not find an acceptable fit of the XMM spectra of the millisecond X-ray pulsar XTE J1751$``$305 using a Comptonization model. While this model seems more appropriate for the spectral results of XTE J0929$``$314, Titarchuk et al. (2002) do not discuss what emission component is responsible for the pulsations.

It is important to realize that both models are based on data with a low-energy bound of 2 keV, which is above the peak flux of the thermal components. The lower energy range of Chandra allows us to make a more robust measurement of $`kT_0`$ for XTE J0929$``$314 than was possible for the RXTE spectrum of SAX J1808.4$``$3658. We conclude that neither model is completely appropriate to explain both our results. More observations of these sources throughout outburst are necessary in order to fully understand their phase-averaged spectral properties. We suggest that in-depth data analysis, including pulse-phase spectroscopy, will provide more insight into the relationship between the various spectral components. The pulse-phase resolved spectral analysis of XTE J0929$``$314 will be presented in a later paper.

We find that the best-fit photon index was significantly different between the Chandra and RXTE fits with $`\mathrm{\Gamma }=1.551.62`$ for the Chandra fits, dependent on the assumed instrumental C abundance, and $`\mathrm{\Gamma }=1.76`$ for the RXTE fits. Since fits with and without the pileup model give consistent parameter values, we reject the possibility that pileup causes the lower Chandra photon index. The combined Chandra and RXTE spectrum of XTE J0929$``$314 is well fit by either a Comptonization $`+`$ blackbody model (which implicitly includes a low-energy turnover from the Comptonization component), or a broken power-law $`+`$ blackbody model, with break energy 1.4–4.4 keV. The high-energy ($`E>3`$ keV) Chandra spectrum is consistent with the RXTE spectral results, indicating that instrumental differences alone do not give rise to the spectral turnover. Based on these results, we suggest that the difference in best-fit photon indices between Chandra and RXTE arises from spectral turnover at low energies. While a power-law model is a reasonable approximation to Comptonization at high energies, extrapolating the power law to low energies leads to a predicted flux that diverges at zero energy. Therefore, the power-law must turnover at some energy for the integrated to flux to be finite.

In addition, the difference in the photon index between the Chandra and RXTE spectral fits is similar to that noted by Miller et al. (2002) between the RXTE and XMM spectral fits of XTE J1751$``$305. From their XMM EPIC spectrum, Miller et al. (2002) found a best-fit power-law photon index of 1.44, while Markwardt et al. (2002) found a best-fit photon index of 1.7–1.9 using RXTE. We interpret the difference between the XMM and RXTE spectral results for XTE J1751$``$305 as evidence of low energy turnover in the broadband spectrum of that source, providing independent support that the turnover is astrophysical in origin.

We would like to thank Chandra Director Harvey Tananbaum and RXTE Director Jean Swank for approving our DDT/TOO observations. In addition, we acknowledge useful discussions with Craig Markwardt, Jon Miller, Mike Nowak, and Dimitrios Psaltis. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. This work was supported in part by NASA under contract NAS8-01129 and grant NAG5-9184.


text:
# The nuance Neutrino Physics Simulation, and the Future

## 1 HISTORY AND OVERVIEW

The unrequited search for proton decay with large underground detectors, begun almost 20 years ago, focused attention on describing atmospheric neutrino reactions in detail. Unlike typical accelerator beams of neutrinos, the spectrum of atmospheric neutrinos covers many decades of baseline and energy. From the lowest-energy contained interactions to the highest-energy upward-going muons, neutrino energies in a broad range ($``$100 MeV to $`>`$1 TeV) are sampled with roughly comparable rates. These features, which would eventually help reveal neutrino mass and oscillation, demand a model of exclusive reactions valid for essentially any interacting neutrino energy.

With the advent of Super-Kamiokande, the statistical precision of atmospheric neutrino measurements increased dramatically and systematic uncertainties (fluxes, cross-sections and calibration) became increasingly important. The same situation is likely to arise in future long-baseline experiments with high-luminosity neutrino beams.

The nuance software package described here is the author’s work-in-progress on the problem of modeling neutrino interactions.111This write-up reflects version 2.000 of the code. Although originally used for atmospheric neutrino interactions, the program was designed to be useful in as many other applications as possible.

## 2 NEUTRINO PHYSICS MODEL

Space does not permit a full elaboration of the cross-section calculations for each reaction channel or graphical comparisons with experimental data, so a schematic outline with references must suffice.

### 2.1 Generalities

Nuance adopts a “divide and conquer” strategy. Models (with varying degrees of sophistication) can be found in the literature for each general class of reaction. By summing the cross-sections and rates of all exclusive channels, and then adding (inclusive) deep-inelastic scattering inside appropriate kinematic limits, the total cross-sections and event rates are obtained.

The units used by the program are MeV, grams (for densities), centimeters and seconds. The user does not directly interact with cross-sections, but these are stored internally in units of $`10^{36}\text{cm}^2`$, or picobarns.

Simulated elementary particles are indexed using PYTHIA/Particle Data Group codes. Application-defined codes are used for nuclei in the few cases they are needed. Particle properties and physical constants are generally taken from PYTHIA internal stuctures by default, although any parameter can be overridden at run-time. Where more recent values are available, they are taken from the Review of Particle Properties. The PYTHIA particle stack is used to build the event in memory and process some decays.

The Fermi gas model is used to simulate the effect of a bound nucleon target (Fermi motion and Pauli blocking). The implementation of this model is slightly different from process to process, however in general bound nucleons are given a uniform initial momentum density up to a (user-specified) maximum value, and a negative binding energy. While different shells (with different momenta and binding energy) are supported, in practice the values $`p_f225\text{MeV/c}`$ and $`E_b=27\text{MeV}`$ (interpolated from a fit to electron-scattering data) are used for all 16 nucleons in the canonical case of $`{}_{}{}^{16}\mathrm{O}`$. A final state nucleon must exceed a threshold momentum value to exit the nucleus and allow the reaction to occur. In the absence of final-state interaction, the (observable) outgoing nucleon momentum must therefore be greater than 225 MeV/c. Discussions at the Workshop made it abundantly clear that this hard cut-off is unphysical, however a more satisfactory presciption remains elusive. In most cases, nucleons below 225 MeV/c will not be detectable, but prudence is advised in interpreting low-energy recoil nucleons.

Lepton masses are never neglected, and the appropriate vector-boson mass is always included in the propagator.

### 2.2 Electron Scattering

The only exactly-calculable (at first-order) neutrino cross-section is unfortunately also the least important for neutrino energies above a few tens of MeV. For completeness, all purely leptonic neutrino and anti-neutrino reactions are treated, including elastic scattering and inverse muon-decay ($`\overline{\nu }_ee^{}\mu ^{}\overline{\nu }_\mu `$). The tree-level cross-section formulæ are given in many particle physics textbooks.

### 2.3 Quasi-elastic Scattering

In nuance quasi-elastic scattering comprises both charged- and neutral-current two-body neutrino reactions with nucleons. The relativistic Fermi gas model of Smith and Moniz provides the general framework for all such processes. To ensure consistency between bound and free targets, free nucleon cross-sections are also calculated using the Smith-Moniz formalism, setting the binding energy equal to zero and expressing the initial nucleon momentum distribution as a delta function at zero.

Identical form factors are used for both free and bound nucleons. The usual dipole parameterization of vector and axial-vector form factors is adopted, with default values $`m_V=0.840\text{GeV}/c^2`$ and $`m_A=1.00\text{GeV}/c^2`$, respectively. For the induced pseudoscalar form factor, a parameterization calculated from lattice QCD (which agrees extremely well with low-$`q^2`$ pion electroproduction data) is used by default, although the simpler pion-pole expression from  can be selected instead. Second-class currents are assumed to vanish, as required by the $`VA`$ electroweak theory and fundamental symmetries.

The total and differential charged-current cross-sections of the Smith-Moniz model agree well with more sophisticated theoretical calculations for neutrino energies above 50-100 MeV where continuum excitation of $`{}_{}{}^{16}\mathrm{O}`$ is dominant.

Neutral-current two-body reactions with nucleons are usually invisible, but have been incorporated into the framework of the Smith-Moniz model and included for completeness. Neutral-current nucleon form factors are specified by the electroweak theory and can be related to those for charged-current reactions.

Similarly, the Smith-Moniz model has been extended to include charged-current, Cabibbo-suppressed hyperon production, following the treatment of Pais to account for the inelasticity of such reactions and the $`|\mathrm{\Delta }I|=\frac{1}{2}`$ rule.

### 2.4 Resonant Processes

For neutrino energies around 1 GeV and above, baryon resonances may be excited and subsequently decay into a nucleon and one or more mesons. The most prominent resonance is the $`\mathrm{\Delta }\text{(1232)}`$, but $`N(1440)`$ and a number of others can also contribute. Rein and Sehgal have shown that these resonance-mediated channels can be described using harmonic oscillator quark wavefunctions, symmetrized under SU(6) (extending the approach of Feynman, Kislinger and Ravndal to account for higher-mass resonant states and the isospin structure of the weak interaction). For hadronic masses above the $`\mathrm{\Delta }\text{(1232)}`$, the data appear to require inclusion of interference between resonances with identical isospin. Nuance adopts the Rein-Sehgal model, modified to account for improved knowledge of the mass spectrum since publication of their original paper, and incorporates non-strange resonances up to 2 GeV. Nucleon form factors for resonance production are assumed identical to those for quasi-elastic scattering, however Rein and Seghal neglect the pseudoscalar form factor of the nucleon, which (in principle) could be important for $`\nu _\tau `$ charged-current reactions.222In practice, deep-inelastic scattering is the dominant $`\nu _\tau `$ charged-current reaction channel for most kinematically interesting neutrino energies.

Rein and Sehgal consider only resonance decays into $`N\pi `$ final states, but additional channels are possible. For completeness, the program attempts to model these by simply adding the additional decay modes so that the total branching fraction for each resonance sums to 100%. The procedure assumes the decay matrix element for other modes is identical to the $`N\pi `$ case, which is probably wrong, but the contribution of these more exotic reactions to the total cross-section is extremely small. Unfortunately, due to the complexity of the Rein-Sehgal calculation, the large number of independent integration variables, and the increased number of channels, the bulk of the program’s computing time is spent on these relatively unimportant reactions.

For reactions on bound nucleons, the initial state is given a uniform Fermi momentum density and a negative binding energy. Reactions in which resonance decay would result in a nucleon below the Fermi sea are Pauli-blocked and do not occur. Kim et al., Singh and collaborators and Marteau et al. have pointed out that the nuclear medium can modify the widths of resonances and allow the final-state reaction $`N^{}NNN`$. In-medium effects on the widths of resonances are not considered by the program, but the “pion-less $`\mathrm{\Delta }`$ decay” reaction (a misnomer, since it seems relevant to $`I=\frac{1}{2}`$ channels as well) can reduce the number of pions produced by 10-50% (the numbers suggested in the literature vary considerably). It is difficult to disentangle $`N^{}NNN`$ from similar processes in which $`N^{}N\pi `$ decay occurs but the pion is subsequently absorbed, and it is not clear to what extent the former reaction is double-counted by final-state interactions. Nevertheless, Super-Kamiokande and K2K data appear to favor some suppression of resonant pion production. Absent a theoretically well-motivated and consistent description of $`N^{}NNN`$ in the Fermi gas model, by default nuance adopts an ad hoc 20% suppression of pion production for $`I_3=\pm \frac{1}{2}`$ reactions and a 10% suppression for $`I_3=\pm \frac{3}{2}`$.333Naïvely, a $`\mathrm{\Delta }^{++}`$ or $`\mathrm{\Delta }^{}`$ can only participate in the $`N^{}NNN`$ reaction with half as many nucleons as $`\mathrm{\Delta }^+`$ and $`\mathrm{\Delta }^0`$. Hopefully interaction between the nuclear- and neutrino-physics communities fostered by the NUINT workshop will help clarify this important effect.

### 2.5 Coherent and Diffractive Reactions

In coherent reactions, neutrinos scatter from an entire nucleus rather than its individual constituents, with negligible energy transfer to the target. Coherent reactions typically produce a forward-going lepton and meson; both charged- and neutral-current reactions are possible, leading to charged or neutral meson production, respectively. As a purely axial-vector process, the cross-sections for neutrinos and anti-neutrinos are equal, and the neutral-current cross-section is half as large as for charged currents.

Nuance implements Rein and Sehgal’s calculation of the cross-section for coherent pion-production in terms of the $`\pi A`$ forward scattering cross-section. Rein and Sehgal’s cross-section agrees with the handful of measurements available, however the lowest-energy data available is at $`E_\nu =2\mathrm{GeV}`$, with an Al target. Using a far more detailed model of the nuclear physics, Kelkar et al. have stressed the importance of nuclear medium effects on the $`\mathrm{\Delta }`$ which mediates these reactions, and predict a dramatic suppression of the coherent cross-section for neutrino energies around 1 GeV. Even with the larger cross-sections predicted by Rein and Sehgal, coherent pion production represents a small fraction of the incoherent, resonance-mediated pion-production cross-section for energies around 1 GeV, so the impact of this effect on the total pion production rate should be 10% or less. Data from the near detectors of K2K, and eventually MiniBooNE, should allow the coherent contribution to be separated from the resonant channels based on the observed single-pion angular distribution.

Diffractive reactions in nuance are analogous to coherent, except free protons rather than $`{}_{}{}^{16}\mathrm{O}`$ nuclei are the targets; the dynamics of the two reactions are identical. Rein’s calculation of the diffractive pion-production cross-section, using a model similar to the coherent case, is implemented in nuance by the same routines which handle coherent reactions. The smaller size of the target reduces the cross-section correspondingly, by approximately a factor $`16^{\frac{2}{3}}`$. Although the nuclear-medium effects cited in  do not apply to free nucleon targets, the cross-section for diffractive single-pion production is dwarfed by incoherent resonant channels and coherent reactions with $`{}_{}{}^{16}\mathrm{O}`$, hence it is included mainly for completeness.

Coherent and diffractive production of vector mesons has been measured by several high-energy experiments, at the few per-mille level compared to the total charged-current rate. Cross-sections for these reactions, which are presently ignored in nuance, have been estimated based on the CVC hypothesis and vector meson dominance and appear to describe the data reasonably well. For completeness, these channels will be added in a future version of the program.

### 2.6 Deep-inelastic Scattering

Deep-inelastic scattering is unique in that it is modeled as an inclusive, rather than exclusive, reaction between a neutrino and the parton consitituents of the nucleon. Elementary formulæ are derived in most particle physics texts, however the usual expressions involve a number of undesirable kinematic approximations, including neglect of lepton (and target, in some cases) masses and the induced pseudoscalar form factor. The unpublished calculation of Roe avoids these limitations, and forms the kinematic basis for the treatment of this reaction in nuance. Additional complications arise at low-energies, where the threshold for production of real heavy quarks (c,b) is important. The usual approach is to impose “slow rescaling”; here the calculation in nuance follows Leader and Predazzi.

Nuance interfaces to the PDFLIB package, allowing one of dozens of nucleon structure function parameterizations to be selected at run-time. By default, the program uses the BEBC option, chosen because it was measured with neutrinos and has a relatively low $`q^2`$ cut-off. Care is required in dealing with structure functions provided by PDFLIB since requests for structure functions with $`|q^2|<|q_{min}^2|`$ result in the values for $`q^2=q_{min}^2`$ being returned. In this case, the behavior of the structure functions as $`q^20`$ is extrapolated according to the vector meson dominance model of PYTHIA.

Shadowing, anti-shadowing, Fermi and “EMC” effects on the structure functions are neglected. Form factors are built from the structure functions assuming the Callen-Gross relation ($`F_2=2xF_1`$) is exact and thereby also neglecting the longitudinal cross-section.

Since the deep-inelastic cross-section includes resonant (and perhaps quasi-elastic) processes with low hadronic mass (which are already included via exclusive channels) the kinematic limits of integration must be chosen to avoid double-counting. As extensively discussed at NUINT’01, in real-life there is no distinct cut-off between “resonant” and deep-inelastic scattering, but rather a smooth transition from scattering off hadrons to scattering off quarks, as the vector boson probe begins to resolve the internal structure of the target. The present version of the nuance attempts to square this circle by integrating the deep-inelastic cross-section over the limits $`|q^2|>1(\mathrm{GeV}/\mathrm{c})^2.\mathrm{OR}.W>2\mathrm{GeV}`$.444The additional constraint $`W>m_N+m_\pi `$ is also imposed, regardless of $`q^2`$. These values are chosen to make the differential cross-section roughly continuous across the artificial boundary at $`W=2\text{GeV}`$, the upper limit for resonant reactions in the program, although it is only partially successful.

Because the deep-inelastic scattering cross-section is inclusive, these channels are different from others in the program, where the final-state particles are uniquely determined by the reaction itself. For deep-inelastic scattering, the interacting quark is determined by a weighting scheme from the LEPTO program and transformed into a different flavor (for charged-current reactions). Parton showers and fragmentation of the outgoing quark are treated by a version of LEPTO which has been extensively modified to ensure conservation of energy and momentum and has kinematic cuts adjusted so $`|q_{min}^2|=10^3(\mathrm{GeV}/\mathrm{c})^2`$ and $`W_{min}^2=2\mathrm{GeV}^2`$. LEPTO itself uses PYTHIA for some of its work.

Because low-mass and small-$`q^2`$ reactions are sometimes generated, the LUND string fragmentation model occasionally fails to produce a final state. In this case KNO scaling is used to determine a final-state using measured multiplicity distributions and fragmentation functions.

Deep-inelastic scattering is undoubtedly the regime where nuance currently lags behind the state-of-the-art in other generators dedicated to high-energy neutrino physics. The NUINT meeting has already borne fruit by bringing the issue of quark/hadron duality to the forefront and eliciting at least one clear prescription for how to remedy these shortcomings. This approach will be implemented in a future version of the program.

### 2.7 Nuclear Processes

The program uses a model of final state interaction in the nucleus originally developed for the IMB experiment. At present the model is specific to $`{}_{}{}^{16}\mathrm{O}`$, although it could be adapted to other nuclei. Primary interactions are assigned a starting position in the nucleus according to the measured density distribution for $`{}_{}{}^{16}\mathrm{O}`$. Hadrons are then tracked through the nucleus in 0.2 fm steps, treating the nucleus as an isoscalar sphere of nuclear matter with radially-dependent density and Fermi momentum. Single-nucleon cross-sections and local density are used to calculate the interaction probability during each step. Interactions resulting in a nucleon with momentum below the Fermi sea are Pauli-blocked and ignored.

Measured cross-sections and angular distributions are used for $`\pi N`$ and $`NN`$ reactions. Angular distributions for elastic reactions are calculated from a global phase-shift analysis of world data. Inelastic reactions involving up to five particles are possible. For three-body $`\pi N`$ reactions, a $`\mathrm{\Delta }I=\frac{1}{2}`$ dominance model is adopted; in other cases, a simple resonance model is assumed. For reactions without available data, cross-sections are inferred from isospin symmetry.

Interactions of kaons with kinetic energy up to 1 GeV are also simulated. Cross-sections and angular distributions for all two-body $`KN`$ reactions are calculated using a partial-wave analysis. Neutral kaons are treated as 50% $`K_S`$ (which decay immediately) and 50% $`K_L`$. $`K_L`$ cross-sections are calculated as an equal mixture of $`K^0\mathrm{and}\overline{K}^0`$, and any interaction has a 50% chance of regenerating a $`K_S`$ which (again) immediately decays.

Hyperon interactions inside the nucleus are not simulated and $`\rho `$ mesons decay rapidly before any interaction is possible. Elastic scattering and pion-production cross-sections compete with decay for $`\omega `$ and $`\eta `$ mesons inside the nucleus.

The IMB cascade model has been tested by running the algorithm in a mode which simulates $`\pi ,p,K{}_{}{}^{16}\mathrm{O}`$ scattering by starting a hadron outside the nucleus and stepping it through. The pion absorption cross-section is tuned to reproduce measurements; for other channels, the results of the simulation reproduce scattering data for hadrons on $`{}_{}{}^{16}\mathrm{O}`$ and $`{}_{}{}^{12}\mathrm{C}`$ extremely well for kinetic energies up to 2 GeV. For pions, the model also agrees with the less detailed (but more elegant) analytical approach of Adler, Nussinov and Paschos.

One flaw in such superposition models was pointed out at NUINT’01, namely the failure to account for “formation zones”, which suppress the final-state interactions of hadrons with energies greater than a few GeV. In the present version of the program, cross-sections are assumed to be constant for energies above 2 GeV.

After ejection of a nucleon recoiling from a neutrino interaction, de-excitation and/or break-up of the residual nucleus can result in emission of one or two few-MeV $`\gamma `$ along with possible evaporation of low-energy nucleons. While the de-excitation products are usually quite low in energy, their signature in a water detector is sufficient to shift the average reconstructed mass of $`\pi ^0`$ produced in neutral current reactions by as much as 5 MeV. The model of de-excitation currently used is specific to $`{}_{}{}^{16}\mathrm{O}`$, and should be generalized.

Once all neutrino interaction products have escaped the nucleus or been absorbed, and the nucleus has de-excited, nuance completes its simulation of the event by calling PYTHIA to process decays of short-lived particles outside the nucleus.

### 2.8 Utilities

Nuance provides additional services for specialized applications, and uses others which are publicly available. The program includes a package for modulating input fluxes under 3-component neutrino oscillations in matter with CP-violation (including a realistic “onion-skin” model of the Earth’s density profile). In practice, one often prefers to generate unoscillated (charged- plus neutral-current) and “fully-oscillated” (charged-current only; $`P_{osc}=1`$, independent of energy) samples, and then reweight them according to different hypotheses after the fact. The program will also generate nucleon decays in a user-specified channel, using the same bound nucleon and final-state interaction model applied to neutrino interactions.

Decays of $`\tau `$ leptons are handled by the TAUOLA package. The polarization of primary $`\tau `$ leptons produced by neutrinos is calculated using the form factors for quasi-elastic and deep-inelastic reactions. In resonant and diffractive processes, the appropriate form factors are unclear and/or difficult to calculate, so the tau is assumed to be completely polarized.

The program can also simulate entering neutrino-induced muons produced in the material outside a detector. For transport of energetic muons to the detector surface, the PROPMU package is used. Note that small-angle Coulomb scattering of muons in flight is neglected for computational efficiency.

A stand-alone utility program (nuplot) is also provided to translate the calculated cross-sections and rates from the internal format of a sequential access ZEBRA/FZ file created by the program into an Ntuple for inspection. Using this utility, energy-dependent cross-sections for any of the many exclusive reactions can be individually plotted.

### 2.9 Validity and Limitations

Unphysical final-state nucleon energies in the Fermi gas model, the $`{}_{}{}^{16}\mathrm{O}`$-specific description of final-state interactions and a number of other problems and simplifications have been discussed above.

Because the program is written in FORTRAN77 and uses ZEBRA for dynamical memory management, single precision calculations are performed in most cases (PYTHIA uses double precision calculations internally). This is an unsatisfactory solution, since it requires considerable care in coding to preserve numerical precision (such problems typically appear for neutrino energies approaching 1 TeV). A future version of the program, using a language with native dynamic allocation (such as FORTRAN90 or C++), will eventually remove this handicap.

## 3 INTERFACE

The program was intended to be flexible and easy to use, although these requirements are frequently at odds, if not mutually exclusive. To allow wider access to the code, it has been posted online at http://nuint.ps.uci.edu/nuance. The code is developed and maintained under the Windows operating system, but is fully compatible with others. The CMT configuration management tool provides a platform-independent mechanism to build the program and its constituent libraries automatically. The main program itself consists of only eight lines of executable code, hence the libraries that do the actual work can be easily slaved to an alternative user-written program.

The program is steered using a ZEBRA/TZ (text) cards file, although frequently-changed options can also be specified on the command line. Default physical parameters are stored in a second cards file named `nuance_defaults.cards`. Values specified in the job-specific file provided by the user override the defaults. The user can also modify any PYTHIA parameter via a data card.

### 3.1 Setup

The user must describe the geometry and composition of the target before running the program. The description syntax is reminiscent of GEANT, but somewhat simpler. A target consists of a one or more volumes (sphere, cylinder or box) which can be nested, but presently must all be centered at the origin. Each volume contains a mixture of one or more materials, with user-specified density. Materials are in turn composed of atoms, which are built from neutrons and protons in shells with specific Fermi momenta and binding energy. Each atom is assumed to include a number of electrons equal to the number of protons.

### 3.2 Fluxes

The program can generate events for a mono-energetic, single-flavor beam of neutrinos, however for most applications the user must also provide a flux of neutrinos. Drivers compatible with several atmospheric flux calculation formats are included.

An accelerator beam can be described via an HBOOK file containing one or more histograms (each associated with a particular neutrino flavor) in units of neutrinos per bin per cm2 per unit luminosity. These histograms must have identical upper and lower limits, but the user may specify whether the limits are interpreted as MeV or GeV. Two scaling factors, one linear and one quadratic, allow the luminosity units and baseline of the beam (respectively) to be changed at run-time. The beam direction in local coordinates is also adjustable. The user may apply multi-quadric smoothing on the input histogram to eliminate binning effects, or directly generate neutrinos from the raw values.

Finally, a simple model of thermal neutrino emission (for supernovæ) is built into program, allowing the user to specify a total luminosity and cooling time for one or more neutrino flavors. In this special case, the fluxes are time-dependent and interaction rates are recalculated as the simulation progresses and the source cools.

### 3.3 Cross-sections and Rates

Having described what is to simulated, the user must calculate neutrino cross-sections and interaction rates. The caculation is time-consuming (a significant fraction of a day on a 2 GHz Pentium machine) due to the multi-dimensional numerical integrations involved, but need only be carried out once for a particular energy range. The cross-sections and rates are stored in a ZEBRA/FZ file for use during event generation; the interaction rates can be recalculated rapidly for a different target size or different flux with identical energy limits by reusing the stored cross-sections. Various classes of reactions, or specific neutrino flavors, may be selectively processed or ignored via flags given on the command line.

Multi-dimensional integration of the rate and cross-section for each channel is performed by the Monte Carlo routine MISER or by cascaded Romberg integration. In both cases, a user-selectable tolerance on the estimated fractional error of the integral (0.3% by default) is checked to determine whether convergence has been achieved. Some kinematic integration variables are transformed internally to accelerate convergence. The integral over neutrino energy necessary to compute total interaction rates can be done linearly or logarithmically. To ensure consistency between calculated rates and simulated data, the same code is used for both cross-section calculations and event generation.

### 3.4 Event Generation

Once cross-sections and rates are calculated and stored, event generation proceeds very rapidly; when generated event vectors are written to disk, the program is I/O-limited. The user indicates the exposure time in units of seconds or years (for atmospheric and supernova fluxes, respectively) or units of luminosity (for neutrino beams). A fixed number of events to generate may alternatively be specified. As in rate calculation, selected classes of reactions or neutrino flavors may be specified or suppressed.

Each reaction channel is tracked separately to determine how much simulation time elapses between events. After determining which channel produces the next event, a neutrino energy is generated from the cumulative differential rate distribution for the reaction stored on the FZ file. A set of transformed kinematic variables is created randomly within the allowed limits and the differential cross-section is calculated. The trial configuration is accepted or rejected by comparing a random number to the ratio of the differential cross-section to the maximum differential cross-section for the same neutrino energy encountered during the rate integration phase (a table of $`\sigma _{max}\mathrm{vs}.E_\nu `$ for each channel is also stored in the FZ file). If a particular set of kinematic variables is rejected, a new set is generated for the same reaction and neutrino energy until an accepted configuration is obtained. Using an accepted set of kinematic variables, four-vectors for the reaction’s primary outgoing particles are generated. These particles are stepped through the target nucleus and/or decayed, if appropriate, and the final outgoing event vectors are written to disk.

For obvious reasons, the program is a voracious consumer of random numbers. The RANLUX package is used everywhere (including external libraries, which are modified to call it instead of their native random number generators). The user should provide a 32-bit integer seed in the steering cards or on the command line to ensure that different jobs produce unique results.

The program writes its results in either (or both) of two formats as requested in the steering cards or on the command line: a human-readable text file or an HBOOK N-tuple. The text file lists not only the initial and final-state tracks, but also the primary outgoing particles prior to nuclear interactions and decays. The reaction channel, elapsed simulation time, neutrino flux at the generated angle/energy and various internal or informational kinematic variables (N-tuple only) are also recorded. Total charged- and neutral-current cross-sections for each neutrino species and target combination are written to the HBOOK file in a separate N-tuple.

## 4 TOWARD A UNIVERSAL MODEL

### 4.1 Motivation

As outlined in the preceding sections, accurately modeling neutrino interactions presents a unique challenge, combining the complexity of the nuclear many-body problem with the formidable uncertainties and sparse experimental guidance typical of neutrino physics. At the same time, discovery of neutrino oscillation has opened a new and inviting window on physics beyond the Standard Model, which ambitious future long-baseline experiments now being planned aim to exploit.

To date, each new experiment has had to solve the problem of simulating neutrino interactions independently. The literature contains a patchwork of disconnected and sometimes incompatible models. Piecing them together involves an equal mixture of intuition, folklore, and guesswork; presently, everyone working on the problem must make their own guesses in a vacuum, largely by trial and error. Without a common point of reference, there is no easy way to test the validity of one’s assumptions and no avenue for one group’s progress or ideas to reach others. In short, there is virtually unlimited room for consolidation and incremental improvement in the state of the art, but no mechanism for it to occur. One can only imagine the chaos which would prevail in $`p\overline{p}`$ or $`e^+e^{}`$ physics if the LUND parton shower and string fragmentation code were not universally available and each experiment were instead forced to recreate something similar on their own, yet this is precisely the present situation in neutrino physics.

Further, as interest in neutrino physics continues to grow, more people enter the field and more future experiments are explored, the need for a general model only increases. On the other hand, the investment of time required to create a reliable simulation is prohibitive for those interested only in exploring future possibilities rather than operating an approved experiment. Most “proprietary” code developed by a particular experimental collaboration tends to be specific to their own running conditions, and many groups are reluctant to release internal software.

### 4.2 The NUINT Working Group

The first NUINT meeting was a watershed in the history of neutrino physics. In addition to bringing together a diverse group from across the spectrum of nuclear and particle physics to offer their specialized insights and expertise, it also united for the first time a large subset of the world’s neutrino simulation experts, many of whom have labored in isolation for years. For me, this opportunity to compare notes and exchange ideas was the high-point of an already-fascinating workshop.

Agreement on the need for more comprehensive and widely available software tools was unanimously expressed, along with a willingness to work together toward them in the future. The first step will be to compare our results in a number of agreed “benchmark” cases. This work will highlight the most important areas of ambiguity and uncertainty, and hopefully motivate renewed theoretical attention on long-neglected but essential points. A central web site (http://nuint.ps.uci.edu) linked to the homepages of each participant will allow the working group and the community as a whole to track the progress of this effort and provide feedback.

For the longer term, the efforts of many will be required to produce a carefully-tested and universal model of neutrino interactions. In addition to purely technical considerations, theoretical guidance and new experimental data will be vital. Still, with the success of NUINT’01 and the promise of renewed and expanded collaboration punctuated and reinforced by future NUINT workshops, it is not too optimistic to hope that within a relatively few years, members of the neutrino physics community will finally have at their disposal a software tool equal to the challenges and opportunities they face.

## 5 Acknowledgements

I would like to thank T. Barszczak, W. Gajewski, A. Habig, D. Kielczewska, P. Lipari, C. McGrew, K. Scholberg, C. Walter, B. Viren and the IMB, Super-Kamiokande and K2K collaborations for important contributions and feedback during the on-going development of the program.


text:
# Solving the riddle of the bright mismatches: hybridization in oligonucleotide arrays

## Abstract

HDONA technology is predicated on two ideas. First, the differential between high-affinity (perfect match, PM) and lower-affinity (mismatch, MM) probes is used to minimize cross-hybridization r2 ; r3 . Second, several short probes along the transcript are combined, introducing redundancy. Both ideas have shown problems in practice: MMs are often brighter than PMs, and it is hard to combine the pairs because their brightness often spans decades r4 ; r5 ; r6 . Previous analysis suggested these problems were sequence-related; publication of the probe sequences has permitted us an in-depth study of this issue. Our results suggest that fluorescently labeling the nucleotides interferes with mRNA binding, causing a catch-22 since, to be detected, the target mRNA must both glow and stick to its probe: without labels it cannot be seen even if bound, while with too many it won’t bind. We show that this conflict causes much of the complexity of HDONA raw data, suggesting that an accurate physical understanding of hybridization by incorporating sequence information is necessary to perfect microarray analysis.

high-density oligonucleotide array (HDONA), mismatch

There are two widespread technologies in use today for performing large-scale mRNA hybridization experiments: spotted arrays and high-density oligonucleotide arrays (HDONAs, a.k.a. GeneChip®) r1 . Such experiments have become popular for assessing global changes in gene expression patterns; they may be used, in a first instance, as screens to identify genes with interesting behaviour on an individual basis; but they also hold the promise to unravel some aspects of the tangled web of transcriptional controls r7 ; r8 . Hybridization array signal is intrinsically "dirty", resulting from compromise to trade quality for quantity, and analysis algorithms therefore need to achieve high levels of noise rejection against the real-world noise observed in the experiments. There is thus a clear need for the early stage algorithms that translate the patterns of light and dark recorded by a laser beam into numbers estimating mRNA concentrations to perform optimally. Any inaccuracies introduced at that level, i.e. loss of signal or false positive assignments cannot be recovered thereafter. In the case of spotted arrays, it seems there is little to do beyond better image analysis; HDONAs however have typically between 20 and 40 probes per transcript, and a function converting those 20-40 numbers into one number has to be supplied r4 ; r9 . As we show below, this task is not trivial, owing to the complex nature of mRNA hybridization and fluorescence detection in this system. HDONA probes are 25-base oligonucleotides grown photolithographically onto a glass surface; about a million different such probes can be synthesized on one chip at current densities. Because 25- mers can exhibit considerable cross-hybridization to a complex background, the system was built on two layers. A "differential signal" approach performs a first level of rejection of spurious signal, by computing the difference between the brightness of a PM probe complimentary to a 25-mer in the target RNA, and a MM probe in which the middle nucleotide has been changed to its complement. From the thermodynamics of DNA-RNA hybrids in solution r10  it was expected that the PM probe should have a higher affinity for the specific target than the MM probe, while cross-hybridization should be roughly equal for both. Second, redundancy was introduced by using several probe pairs corresponding to distinct 25-mers along the length of the transcript (see Figure 1).

But these ideas do not translate that easily from hybridization in solution to HDONAs. An issue long noticed was the large number of probe pairs for which the single mismatch brightness was higher than the perfect match up to a third of all probe pairs in some chip models 6. This was easy to notice since early versions of the default analysis software would not take this matter into account, and therefore some gene concentrations were reported as negative. Why this would happen has been the cause of much speculation. A two-dimensional plot of PMs vs. their MMs shows that their joint probability distribution appears to have two branches, and it was suggested that sequence specific effects are playing a crucial role 6. But in the absence of sequence information for the probe pairs, this couldn’t be verified. However, Affymetrix has recently released the necessary data for addressing the problem explicitly.

We show in Figure 2 joint probability distributions of PMs and MMs, obtained by taking every probe pair in a large set of experiments, and binning them to obtain two- dimensional histograms. We did this twice, computing two separate probability distributions which we then superimposed: in red, the distribution for all probe pairs whose 13th letter is a purine, and in cyan those whose 13th letter is a pyrimidine. The plot clearly shows two very distinct branches in two colors, which correspond to the basic physical distinction between the shapes of the bases: purines are large, double ringed nucleotides while pyrimidines are smaller single ringed ones. This underscores that by replacing the middle letter of the PM to its complementary base, the situation on the MM probe is that the middle letter always faces itself, leading to two quite distinct outcomes according to the size of the nucleotide. If the letter is a purine, there is no room within an undistorted backbone for two large bases, so this mismatch distorts the geometry of the double helix, incurring a large steric and stacking cost. But if the letter is a pyrimidine, there is room to spare, and the bases just dangle. The only energy lost is that of the hydrogen bonds.

So the existence of two branches agrees with basic hybridization physics, but it still does not explain why the MMs are actually brighter than the PMs in many sequences with a purine middle letter. To understand this we will perform a finer level of analysis, concentrating momentarily only on the PM sequences. It has been pointed out that the PMs within a probeset are very broadly distributed, typically spanning two decades or more. We can try to observe whether this breadth is similarly sequence- dependent, by fitting the brightness B of PM probes (divided by the estimated RNA concentration \[RNA\]) against their own sequence composition:

$$\mathrm{log}\left(B/[RNA]\right)=\underset{sp}{}L_{sp}A_{sp}$$

where $`s`$ is the letter index (ACGT) and $`p`$ the position (1-25) on the 25-mer; $`L`$ is a Boolean variable equal to 1 if the symbol $`p`$ equals $`s`$, and thus $`A`$ is a per-site affinity. More accurate models would include stacking energies by looking at consecutive letters (bonds); while this contribution is important for hybridization experiments in solution r11 ; r12 , we found that it does not improve the fit substantially. On the other hand, we were surprised to discover that the key improvement comes from introducing position dependent affinities, as opposed to affinities that would depend only on the total number of occurrences of each letter. The fitted per-site affinities are shown in Figure 3, note the strength of letter specific contributions: changing an A to a C in the middle of the sequence would change the brightness of the probe by 250on mouse, drosophila, and yeast arrays lead to virtually identical affinity curves as those shown in Figure 3. Besides providing insight into physical aspects of hybridization, the fitted affinities bear an important practical value as they permit to effectively reduce the breadth of the probeset brightnesses, therefore improving the signal-to-noise ratio of probeset averages (used for instance in absolute concentration estimates). In numbers, the variance in 96by the fit is subtracted, and the reduction is larger than a factor of 2 for 65probesets. An interesting aspect of the above fits is the asymmetry of A vs. T (and G vs. C) affinities, which is shown more clearly in Figure 4.

The obvious culprits for this effect are the labels, namely, the standard protocol recommended by Affymetrix entails labeling the cRNA with biotinilated nucleotides more specifically, U and C, the pyrimidines. This suggests a rather simple explanation, namely, that the biotinilated bases somehow impede the binding; the effect diminishing to zero toward the probe edges, where the double strand breathes enough to be able to accommodate the linkers, and being maximal near the center, where the largest disruption would be effected. This would cause a catch-22 in terms of obtaining the maximal fluorescence: if a sequence has too few bases that can be labeled, it will not shine, even if it binds strongly, while if it has too many labels it will not shine because it does not bind. But this catch-22 has a curious loophole: the optimal region to have the fluorophores should then be outside the 25-mer: since the cRNA being hybridized is usually longer. Figure 4 confirms this: when including the contribution to brightness from sequence composition outside the 25-mer we find the pyrimidine contribution to be strictly positive.

Interference with binding by the biotinilated bases also solves the MM\>PM riddle. As we saw before, a purine in the middle of the PM implies a gap between the two nucleotides on the MM probe; thus one could conjecture that this gap permits the linker between nucleotide and biotin not to interfere with the binding. This conjecture is quantitatively compatible with the data: according to Figure 4, the energetic penalty for a pyrimidine in the middle of the sequence is 0.2 in log10 units (about 0.5 kBT), which is comparable to (and bigger than) the excess brightness of the MMs in the purine (red) lobe of Figure 2. Indeed, the median excess brightness of the MM for the red probes is 0.1 in log10 units. In other words, when considering the effective contribution of a middle bond to brightness, a G-C\* bond on the PM probe is dimmer than a C-C\* bond on the MM, which in turn is dimmer than a C-G bond on the PM. Here * denotes a labeled nucleotide on the mRNA strand.

A microarray experiment carried out for a biological study provides nonetheless a quarter of a million measurements in hybridization physics. This information may be used to probe and understand the physics of the device, and indeed if an accurate enough picture emerges, it shall lead to substantial improvements in data quality. We have shown how the basic physics of the detection process in HDONAs percolates into the statistics, resulting in statistical anomalies affecting the data thereafter, and which need to be taken into account in order to optimize the experiments. Microarrays are one out of many high-throughput techniques being developed and brought to bear in important problems in Biology today. While it is usually emphasized that they pose similar analytical challenges in terms of pattern discovery, mining and visualization, our work exemplifies that in order to reach a level where analysis can be abstracted to such heights, one should be positive to understand in some detail the physics of the instrument and how it affects the data. We thank Herman Wijnen, Edward Yang, Nila Patil, Coleen Hacker and Adam Claridge-Chang for helpful discussions. Current address (MM): The Abdus Salam International Centre for Theoretical Physics, Strada Costiera 11, Trieste I-34100.


text:
# Hydrogen-induced Unzipping of Single-Walled Carbon Nanotubes

## Abstract

The chemisorption of atomic hydrogen on the single-walled armchair carbon nanotube is studied with ab initio calculations. A single H atom is found to be chemisorbed on both the inside and outside wall of the nanotube. The binding energy of H adsorption at the exterior of the nanotube is much greater than that at the interior of the nanotube. For the first time, we predict that two rows of H atoms chemisorbed on selective sites exterior to the nanotube can break the nearest-neighbor C-C bond of the nanotube through the concerted formation of C-H bonds, leading to the unzipping of the nanotube wall. We provide insights into the underlying electronic structure responsible for the H-induced unzipping of the nanotube, lending strong support to the recent experimental observations for the coalescence of single-walled nanotubes in the presence of atomic hydrogen. Interestingly, H atoms chemisorbed inside the nanotube do not lead to the breaking of the C-C bonds.

Carbon nanotubes have been attracting considerable interest due to their unique electronic and mechanical properties since their discovery about a decade ago dresselhaus1 . In particular, there have been intense studies in evaluating the capability of carbon nanotubes as hydrogen storage material for clean energy sources chen ; liu . It was found that hydrogen could be stored in bundles of the single-walled nanotubes up to 5-10 wt.% at pressures less than 1 bar near room temperature dillon . The mechanism, for the most part, is attributed to physisorption of H2 molecules both inside the tubes and within the interstitial regions among the nanotubes at room or low temperatures dillon ; ye . However, recently Liu et al. observed that, after treating carbon nanotubes with H2 gas under high pressure, there was residual H2 during the desorption cycle liu , suggesting the presence of chemisorption in the process. On the other hand, chemisorption can also take place during electrochemical storage processes, in which H2 molecules are broken into H atoms with the application of metal catalysts or electrochemical techniques haufler ; becker . On the theoretical side, first-principles calculations have confirmed the dissociative chemisorption of H2 on two adjacent nanotubes in solid phase, proceeding with the breaking of the H-H bond concerted with the formation of two C-H bonds on two adjacent nanotubes chan . The chemisorption of atomic H on single-walled nanotubes has also been observed from self-consistent-charge-density tight-binding calculations lee2 ; lees .

Recently, Nikolaev and coworkers discovered a remarkable phenomenonsmalley : under atomic H atmosphere, single-walled armchair nanotubes annealed up to 1500C coalesce with neighboring tubes, resulting in larger nanotubes with twice and occasionally three times the diameter of the original ones. Based on these observations, the authors proposed a H-activated coalescence mechanism, in which the gas-phase H atoms attack the side of neighboring nanotubes, breaking the C-C bonds, and producing defective sites on the adjacent nanotubes. Once these adjacent defects are formed, the strong thermodynamic force, resulting from the released strain energy in forming larger tubes, drives the two neighboring smaller tubes to join together. Coalescence of clean single-walled nanotubes has also been observed under electron irradiation conditions at elevated temperatures terrones .

In this paper, we present a systematic ab initio study of chemisorption of H atoms on a single-walled (6,6) armchair carbon nanotube. For the first time, we provide theoretical evidence for the H-induced unzipping of the nanotube. Our theoretical calculations lend strong support to the recent experimental observations of H-activated coalescence of armchair nanotubes smalley , revealing the electronic mechanism responsible for this remarkable phenomenon.

The ab initio calculations we performed are based on density-functional theory with the CASTEP implementation and ultrasoft pseudopotentials vanderbilt . The Perdew-Burke-Ernzerhof gradient-corrected functional perdew  was used for the exchange-correlation potential. The energy cutoff for the plane-wave basis was set to be 300 eV, yielding a convergence for the total energy better than 1 meV/atom. For the reciprocal-space integration we have used eighteen $`k`$-points along the direction corresponding to the nanotube axis using the Monkhorst-Pack scheme monkhorst . Results have been obtained for the fully relaxed geometries including all atoms and the lattice constant of the supercell along the tube axis.

The single-walled (6,6) armchair carbon nanotube is modeled with a supercell of dimension 20Å$`\times `$ 20Å$`\times `$ 2.44Å. This corresponds to two layers of C atoms perpendicular to the tube axis with twelve atoms per layer. We have studied three different H coverage, i.e., 1, 2 and 24 H atoms per unit cell. We should point out that due to the periodic boundary conditions along the tube axis, 1, 2 and 24 rows of H atoms are actually being simulated in the calculations. The diameter of the fully relaxed pure nanotube is found to be 8.17 Å. The binding energy per H atom,

$$E_b=1/n[E_{tot}(tube)+E_{tot}(nH)E_{tot}(tube+nH)]$$
(1)
is calculated in terms of the total energy of the pure nanotube $`E_{tot}(tube)`$, the total energy of the $`n`$ (=1, 2, 24) H atoms $`E_{tot}(nH)`$, and the total energy of the nanotube with $`n`$H atoms, $`E_{tot}(tube+nH)`$. Since $`E_{tot}(nH)`$ is calculated using the same supercell geometry as in the ($`tube+nH`$) system, the spurious adatom-adatom interactions along the tube axis are subtracted. In Table I we list values of the binding energy per H atom at the interior and exterior of the nanotube for one H atom, a pair of H atoms on the same and adjacent layers perpendicular to the nanotube axis, as well as the fully hydrogenated case. The positive value of $`E_b`$ in all cases indicates that the chemisorption is exothermic and hence energetically stable.

First, we consider the H coverage of one atom per unit cell, with the H atom adsorbed either at the interior or exterior of the tube. By placing the H atom at various initial positions away from the nanotube wall, we find that there is no energy barrier for H chemisorption at both inside and outside of the tube. Fig. 1(a) shows the relaxed atomic structure with the H atom adsorbed at the exterior to the tube. We label the nearest-neighbor C atom to the H by C1, and the nearest-neighbor C atom to C1 at the same layer by C2. The calculated C1-H bond length (1.11 Å) is the same for the H atom at either side of the nanotube wall, and it is close to the corresponding value of 1.10 Å in a CH4 molecule. However, there is a significant difference in the binding energy for H at the interior (0.49 eV) and the exterior (1.77 eV) of the nanotube wall. This is due to the fact that the chemisorption of H results in a transition from $`sp^2`$-like bonding in the pure nanotube to $`sp^3`$-like bonding in the presence of H; and the latter cannot be fully formed upon adsorption to the inner wall because of the unfavorable bonding angle tada . The Mulliken analysis of bond population mulliken ; sanchez ; segall  shows that 0.33 electrons is transferred from the H atom to its nearest-neighbor C1, in agreement with previous calculations ciraci  for the fully hydrogenated nanotubes. The C1\- $`s`$ and $`p`$ orbitals acquire 0.15 and 0.18 electrons, respectively. It is important to note that for the pure nanotube the overlap population between the nearest-neighbor C atoms on the same and two adjacent layers is 0.87 and 2.22, respectively, indicating a stronger bonding in the latter case sanchez . This large difference in overlap population between the two types of nearest-neighbor C bonds will in turn affect the propensity of H-induced breaking of the C-C bonds. The chemisorbed H results in an increase of the C1-C2 bond length from 1.42 Å (in the pure nanotube) to 1.53 Å, with the latter value being typical for $`sp^3`$ C-C bonds. As a result, the C1-C2 overlap population decreases from 0.87 to 0.63, indicating that the C-C bond is weakened. In order to gain more insights into the change of the electronic structure accompanying the H chemisorption process, we have calculated the angular momentum- and site- projected density of states (DOS) for some representative cases (Fig. 2). For the pure nanotube, we have confirmed the metallic behavior of the (6,6) armchair nanotube with a low density of states at the Fermi energy ($`E_F`$) dresselhaus1 . From Fig. 2(a), we find that the H chemisorption induces the opening of a band gap ($``$ 5 eV) for the C1 atom, which is the result of the large charge transfer from the H atom to C1, filling up the C1 $`p`$ valence band. On the other hand, it is interesting to note that the C2 atom exhibits a large narrow peak across $`E_F`$ which overlaps with the bonding states of the H atom. This peak, reminiscent of a Van Hove singularity, arises from the degeneracy of the two linear bands crossing $`E_F`$ in the pure nanotube, which is in turn lifted under the local distortion induced by H, thus forming two sub-bands whose extremes fall at EF. Fig. 3(a) shows the valence charge density contour plot for a plane containing the C atoms on the same layer. One can see the directional $`sp^2`$-like covalent bonding between the C atoms, the $`sp^3`$-like bonding between C1 and H, and the weakening of the C1-C2 $`sp^2`$ bond.

The result of the H-induced weakening of the C-C bond invites an interesting question: Could the chemisorption of another H atom on the nearby C atom lead to the formation of two H-C bonds which will in turn further decrease or even break the C1-C2 bond, resulting in the unzipping of the entire nanotube? Thus, next we consider the case of a pair of H atoms chemisorbed on the C1 and C2 atoms in the exterior and interior to the tube. We have examined two different adsorption sites for the second H atom, either at the same or the adjacent C layer. First, we present results with the interatomic axis between the H atoms parallel to the C1-C2 bond. The binding energy per H for this case is 0.37 eV (inside the tube) and 3.0 eV (outside the tube), respectively. As in the single H case, the pair of H atoms prefers to be chemisorbed at the exterior of the nanotube. Interestingly, the binding energy per H atom for the pair outside the nanotube has dramatically increased compared the corresponding value for the single H case, while the charge transfer from each H to its nearest-neighbor C remains about the same (0.30 electrons) as in the one H case. This charge transfer leads to an ionic-like repulsion between the H atoms and between C1 and C2. Figure 1(b) shows the relaxed atomic structure for the case of two H atoms adsorbed at the exterior of the tube on the same layer. The most remarkable feature in this figure is the breaking of the nearest-neighbor C1-C2 bond. The equilibrium H-C, C1-C2, and H-H bond lengths are 1.09 Å, 4.07 Å, and 2.94 Å, respectively. Due to the periodic boundary conditions along the tube axis, the result implies scission of the nanotube wall along the tube axis, i.e., the unzipping effect necessary for the experimentally observed coalescence phenomenon. The overlap population of C2-H is 0.71, which is as large as that for the nearest-neighbor C-C bond at the same layer in the pure nanotube. Fig. 2(b) shows the DOS for one of the hydrogenated C atoms, C2 (see also Fig. 1b) and its nearest-neighbor, C3. The DOS for C2 is similar to that of C1 in Fig. 2(a). Note the narrowing of the H $`s`$ band which exhibits bonding and antibonding states separated by an energy gap of about 5 eV, indicates that the H atom is stabilized due to the hybridization between the $`s`$ states of H and the $`p`$ states from the nearby C2 atom. From Fig. 3(b), one can see the $`sp^3`$-like bonding between the C and the H atoms, and the breaking of the C1-C2 bond.

We have also studied the case of two H atoms chemisorbed on the two nearest-neighbor C atoms located at adjacent layers shown in Fig. 1(c). The binding energy per H atom is found to be 0.73 eV in the interior and 2.52 eV in the exterior of the tube, respectively. In contrast to the previous case, we find that the nearest-neighbor C-C bond does not break. The equilibrium H-C, C2-C3, and H-H bond lengths are 1.10 Å, 1.50 Å, and 2.29 Å, respectively. There is a similar amount of charge transfer from the H atoms to each nearest-neighbor C atoms with an overlap population of 0.74. However, in contrast to the previous case, the C2-C3 overlap population of 1.59 is reduced considerably relative to that of 2.2 in the pure nanotube. By employing a larger supercell (four layers per cell), we find essentially the same results, i.e., the C-C bond remains unbroken in the presence of H atoms. This result can be understood from the fact that the nearest-neighbor C-C bonds on adjacent layers are stronger than those on the same layer, as alluded earlier.

Interestingly, the C-C bond can not break when the two H atoms chemisorb inside the nanotube regardless of the adsorption site. The C-C bond deforms nevertheless, with the bond length increasing from 1.42 Å (pure tube) to 1.52 Å. This again originates from the reduced propensity for $`sp^3`$ bonding at the interior of the nanotube tada . Overall, the internal surface of the nanotube is found to be less reactive than the external surface, in agreement with the general consensus haddon ; cho . Finally, our calculations for the fully endo- and exo-hydrogenated nanotube find no breaking of the C-C bonds and a lower binding energy per H atom of 0.30 eV and 0.84 eV, respectively. We have also carried out calculations for a H2 molecule approaching the nanotube wall, both from the exterior and the interior, with the H-H bond perpendicular to the nanotube wall. In agreement with previous calculations chan , we find that no chemisorption occurs and the H2 molecule remains intact.

In conclusion, we have studied the chemisorption properties of atomic H on a single-walled (6,6) armchair nanotube. We find that H atoms bind strongly to the nanotube through $`sp^3`$ bonding, with the binding energy at the exterior of the tube being much greater than that at the interior. For the first time, we predict that a pair of H atoms chemisorbed on two nearest-neighbor C atoms on the same layer outside the nanotube catalyzes the breaking of the C-C bond and leads to the unzipping of the nanotube. This H-induced unzipping mechanism lends strong support to the recent experimental observations for the coalescence of single-walled nanotubes in the presence of atomic H. On the other hand, the C-C bond does not break when H atoms are adsorbed at the interior of the tube.

We acknowledge support from the U. S. Army under Grant Nos. DAAD19-00-1-0049 and DAAG55-97-1-0247.