text:
# BeppoSAX observation of PSR B1937+21

## 1 Introduction

The various X-ray satellites, from ROSAT to CHANDRA, observed about a dozen millisecond pulsars (MSPs) (see e.g. Becker 2001, Becker & Aschenbach 2002). Such observations demonstrate that the X-ray emission from MSPs is mainly not of thermal origin and those with the hardest spectra appear to be objects with strong magnetic fields $`B_L`$ at the light-cylinder radius $`R_L`$ ($`R_L=cP/2\pi `$ where $`P`$ is the spin period of the MSP; see Saito et al. 1997, Kuiper et al. 1998, Kuiper et al. 2000, Becker 2001). However results, in the soft “ROSAT energy band” 0.1–2.4 keV from a Chandra observation of 47 Tuc, suggest that significant thermal emission (at energies $`kT0.2`$ keV) can be produced from MSPs (see Grindlay et al. 2002 for details). This appears not to be the case for PSR B1821$``$24 (see e.g. Becker et al. 2003). Correlation between spin-down energy loss and X-ray luminosity was investigated by several authors (e.g. Verbunt et al. 1996, Becker & Trümper 1997, Takahashi et al. 2001, Possenti et al. 2002), suggesting that the claimed law $`L_X`$(2–10 keV) $`\dot{E}^\gamma `$, with $`\gamma `$ in the range 1.0–1.5, is valid for MSPs and “ordinary” pulsars in the same way.

Unfortunately, good spectral and temporal informations exist only for about half of the targeted sources. In fact (with the exception of the data collected by Rossi-XTE and, partially, by ASCA and BeppoSAX) MSP X-ray observations are usually affected not only by low statistics but sometimes also by insufficient time accuracy to perform detailed periodicity and timing analyses.

PSR B1937+21, with a period $`P1.56`$ ms, is the first and still the fastest rotating MSP known. Its radio pulse profile is double peaked with phase separation of about 0.52. In spite of its low surface magnetic field strength of $`B_S=3.2\times 10^{19}(P\dot{P})^{1/2}=4.1\times 10^8`$ G, its magnetic field at the light-cylinder is the highest of all known pulsars: $`B_L=B_S(R_S/R_L)^3=2.97\times 10^8\dot{P}^{1/2}P^{5/2}1\times 10^6`$ G (it is $`R_L7.4\times 10^6`$ and the neutron star radius $`R_S=10^6`$ cm), very similar to that of the Crab pulsar. With a spin-down luminosity of $`\dot{E}=4\pi ^2I\dot{P}/P^31\times 10^{36}`$ erg s-1 ($`I=10^{45}`$ g cm2, momentum of inertia), the derived spin-down flux density $`\dot{E}/(4\pi (d/3.6\mathrm{kpc})^2)7\times 10^{10}`$ erg cm-2 s-1 is similar to that of PSR B1821$``$24 but it is 10 times larger than that of PSR J0218+4232, i.e. the other two MSPs detected in X-rays with good statistics which have high luminosities and show non-thermal emission only.

Analysis of a sensitive observation performed by Rossi-XTE has revealed that the X-ray peaks (see below) are almost perfectly aligned with the radio giant pulses (Cusumano et al. 2003). Such giant pulses are short bursts of emission with a flux density exceeding that of the whole profile by a factor of 10 to 100 or more. Having been discovered first for the Crab pulsar (Staelin & Reifenstein 1968), they have also been detected for PSR B1937+21 (Wolszczan et al. 1984, Backer 1995), and more recently for PSR B1821$``$24 (Romani & Johnston 2001), PSR B0540$``$69 (Johnston & Romani 2003) and PSR B1112+50 (Ershov & Kuzmin 2003). The observations suggest that the radio giant pulses are related to the high energy emission of these pulsars, possibly originating in outer gaps (Romani & Johnston 2001). Interestingly, Kinkhabwala & Thorsett (2000) report that the giant pulses of PSR B1937+21 appear in windows of $`55÷70`$ $`\mu `$s located after the main and interpulse peaks, with a typical duration of less than 10 $`\mu `$s. Thereby, the location of the X-ray peak coincides with the giant pulses and not with the normal radio profile, strongly supporting the idea about a common origin of giant pulses and high energy emission.

In this paper we present the results of the temporal and spectral analysis of a BeppoSAX observation of PSR B1937+21. Compared to those performed with ASCA (Takahashi et al. 2001) our observation exploits the better instrumental sensitivity and a longer exposure time. This allows the detection of the pulse profile interpulse and source photons down to energies of 0.5 keV and enables us to put constraints on the $`N_\mathrm{H}`$ toward the pulsar. Our results are compared with those obtained with ASCA and, more recently, with Rossi-XTE (Cusumano et al. 2003).

## 2 Observation and spatial analysis

BeppoSAX observed PSR B1937+21 on May $`1^{\mathrm{st}}`$ 2001. Pointing coordinates were those from the radio timing. The $`200`$ ks observation time resulted in a total effective on-source time of $`78.5`$ ks and $`34`$ ks with the MECS (1.3–10 keV; Boella et al. 1997) and LECS (0.1–10 keV; Parmar et al. 1997), respectively. To extract the source photons from LECS and MECS we followed two alternative ways:

* extraction of the photons from the circular region and in the energy range which optimize the signal to noise ratio: they were $`3^{}`$, 0.5–8 keV for the LECS and $`2^{}`$, 1.3–10 keV for the MECS, resulting in 106 and 385 collected photons, respectively. The local background (in a circular region in the field of view $`10^{}`$ away from the pulsar) was compared to that obtained using archival data of blank sky observations. We found that the local background is $`12\%`$ higher for both LECS and MECS; background photons were 50 and 123, respectively.
* a Maximum Likelihood (ML) approach to extract the number of counts assigned to the pulsar taking into account the presence of other sources and the background simultaneously (see e.g. Kuiper et al. 1998).

## 3 Temporal analysis

The arrival times of the extracted photons were converted to the Solar System Barycentric (SSB) frame using the BARYCONV111see www.asdc.asi.it/bepposax/software/saxdas/ code and then searched for periodicities by folding the MECS data with trial frequencies around the radio one222 Folding of the lower statistics LECS data did not reveal any modulation at $`>2\sigma `$ level.. The timing ephemeris listed in Table 1 was obtained from observations with the Effelsberg 100-m radio telescope and archival Arecibo data. Details are reported in Cusumano et al. (2003). We found a clear signal, but the corresponding frequency $`f_X=641.92825315`$ Hz deviated from the radio value by $`(f_Rf_X)7.6\times 10^6`$ Hz. This difference is about ten times the statistical error ($`8\times 10^7`$ Hz). It could be ascribed to incorrect photon time markers, likely due to long term degradation of the on-board clock. This effect is relevant for long duration observations only (see Nicastro et al. 2002 for more details).

Figure 1 shows the folded light curve of the MECS data together with the radio one obtained with the 100-m Effelsberg radio-telescope. The X-ray profile is double peaked: the primary peak (P1) is detected at $`10\sigma `$ level, the secondary peak (P2) at $`4\sigma `$. We performed a Gaussian fit on the two peaks. The derived phase widths are $`\sigma _{\mathrm{P1}}=0.030`$ and $`\sigma _{\mathrm{P2}}=0.024`$, respectively. This values translate in $`100`$ $`\mu `$s FWHM, fully compatible with the instrument time resolution limit. The phase separation of the two peaks (primary to secondary) is $`0.52\pm 0.02`$ (90% confidence), where the error is derived from the statistical uncertainties in the peak position gaussian fits: $`\pm 0.007`$ and $`\pm 0.010`$ for P1 and P2, respectively, Figure 2 shows the X-ray pulse profile in 3 energy ranges. Non-detection of the interpulse in the ASCA data is explained by their lower statistical significance (227 total photons in 1.7–6.5 keV).

To quantify the pulsar DC component, we applied the bootstrap method proposed by Swanepoel et al. (1996), which allows us to estimate the off-pulse level (sky background plus unpulsed source component) in the pulse profile. We find that over the total energy range 1.3–10 keV the source has a pulsed fraction333It is $`f=`$ the percentage of total counts above DC level with respect to the net total counts in the full profile above the background level (determined in the spatial analysis). of $`f=85\pm 5\%`$ reducing to $`54\pm 7\%`$ at energies 1.3–4 keV. In the 4–10 keV band the detected photons are $`100\%`$ pulsed. Changing the folding initial phase causes the pulse profile to change slightly but the detection level of the two peaks and the pulsed fraction remain well within the statistical error. The ASCA pulsed fraction of 59% is compatible with our results.

There is an indication for an increase of the strength of P2 compared to P1 with energy. Using the bins marked in Fig. 1 as on-pulse phases, $`P1/P2=5.13\pm 1.69`$ for the band 1.3–4 keV and $`P1/P2=2.99\pm 0.84`$ for the band 4–10 keV.

A joint analysis of radio data (Jodrell Bank) and X-ray data (ASCA) led Takahashi et al. (2001) to claim an alignment of the main X-ray peak with the radio interpulse. The BeppoSAX on-board clock absolute accuracy does not allow us to align our profile to the radio one. However we noticed earlier (Nicastro et al. 2002) that the radio and X-ray peaks have both a phase separation of $`0.52`$ (primary to secondary), then leading us to the conclusion that their result is likely incorrect. Results from a Rossi-XTE observation has confirmed this doubt showing that the X-ray peaks are almost perfectly aligned with the radio giant pulses (see Cusumano et al. 2003).

## 4 Spectral analysis

In spite of the lower exposure time and sensitivity of LECS we were able to detect a significant “unpulsed” signal in the 0.5–8 keV range for this instrument. Inclusion of these softer photons in the spectral fit allowed us to better constrain the $`N_\mathrm{H}`$ toward the pulsar. In grouping the LECS and MECS data we applied the condition to have at least 20 and 30 counts per bin, respectively. The best fitting model (in the detection range 0.5–10 keV) was an absorbed power-law. A Black Body model did not fit the total flux data satisfactorily ($`kT1.2`$ keV and $`\chi _n^2=1.5`$, 13 d.o.f.). We then tried to fit the data with a broken power-law and a power-law + Black Body. For the latter case we also tried fixing the value of kT to various values in the range $`0.1÷0.5`$ keV as higher values are unlikely for an MSP (Zavlin et al. (2002) found $`kT0.2`$ keV for PSR J0437$``$4715 and Grindlay et al. (2002) report $`kT0.22`$ keV for the MSPs in 47 Tuc). The F-test gives chance improvement probabilities of 90% for the broken power-law model and values around 50% for the power-law + Black Body, then we can reject these models. In the energy bands 0.1–2.4 and 2–10 keV we can estimate *absorbed* flux upper limits of $`1\times 10^{14}`$ and $`5\times 10^{14}`$ erg cm-2 s-1 for a Black Body component with $`kT=0.5`$ keV. On-pulse analysis, subtracting the off-pulse counts, was also performed; phase bins used to extract the pulsed photons are marked in Fig. 1. Fit of the 1.3–10 keV MECS data alone gave $`N_\mathrm{H}=(1.2\pm 0.9)\times 10^{22}`$ cm-2 and $`\alpha =1.7\pm 0.4`$ ($`\chi _n^2=1.06`$).

In Table 2 the best fit results (for the absorbed power-law model) are reported; in Fig. 3 the total (LECS+MECS) and pulsed spectra are shown together with the unabsorbed power-law and the fit residuals for both the total and the pulsed spectra. We also fitted the photons collected in the off-pulse region (DC level) which gave, though at lower statistics, a consistent $`N_\mathrm{H}`$ and a photon index of $`3.3\pm 0.5`$. We note that the spectral index of 1.2 of the pulsed emission is harder than that obtained from the total emission. This is a typical value for magnetospheric emission and is fully consistent with the value obtained by Cusumano et al. (2003) using Rossi-XTE data. The total unabsorbed and absorbed fluxes of PSR B1937+21 in the 2–10 keV range are $`4.1\times 10^{13}`$ and $`3.5\times 10^{13}`$ erg cm-2 s-1, respectively. The derived (unabsorbed) luminosity is $`L_X=5.0\times 10^{31}\mathrm{\Theta }`$ ($`d/3.6`$ kpc)2 erg s-1 while the X-ray conversion efficiency is $`\eta =L_X/\dot{E}=4.5\times 10^5\mathrm{\Theta }`$, where $`\mathrm{\Theta }`$ is the solid angle spanned by the emission beam. This result is confirmed by the Rossi-XTE 2.5–17 keV spectrum (Cusumano et al. 2003). We also note that the fit parameters are consistent with the non-detection (Verbunt et al. 1996) in the ROSAT HRI energy range (0.1–2.4 keV).

## 5 Discussion

### 5.1 Neutral hydrogen column density

Galactic $`N_\mathrm{H}`$ can be estimated from radio/IR measurements in three different ways: 1. from the DM by assuming 10 H atoms for each e-; 2. from H i measurements, assuming there is no hydrogen at molecular state; 3. from dust extinction using the relation $`N_\mathrm{H}=1.79\times 10^{21}A_V`$ of Predehl & Schmitt (1995). When there are no direct pulsar measurements or three dimensional distribution models, points 2. and 3. make use of models of the integral content of H along a particular line of sight. They then only give upper limits for the atomic and total hydrogen content, respectively.

For PSR B1821$``$24 placed in the globular cluster M28, at a distance of $`5.5`$ kpc and galactic latitude of $`5.^{}6`$, one finds $`N_\mathrm{H}`$ $`1.6\times 10^{21}`$ cm-2 which is consistent with the estimates from Galactic H i, dust and pulsar dispersion measure (Becker et al. 2003). Similarly consistent results are obtained for PSR J0218+4232. However, because of variable e-/$`N_\mathrm{H}`$ ratios in the various Galactic directions and the presence of more o less dense gas clouds along the line of sight, discrepancies would not be surprising.

The column density we derive from the power-law fit is 10 times greater than the radio DM derived one and is only marginally consistent (within $`2\sigma `$) with the H i derived Galactic value $`1.2\times 10^{22}`$ cm-2 (Dickey & Lockman 1990). In addition, it is apparently totally inconsistent with the $`7.4\times 10^{22}`$ cm-2 obtained from the Schlegel, Finkbeiner & Davis (1998; SFD) dust extinction maps (Schlegel et al. (1998).

Being the line of sight to PSR B1937+21 tangent to the Galactic spiral arms ($`l57.51`$, $`b0.29`$, see e.g. Fig. 3 of Watson et al. 2003), one can invoke a particularly low e-/$`N_\mathrm{H}`$ ratio in that direction due to the lack of ionizing sources. H i maps report the integral Galactic column density but do not account for hydrogen in molecular form and then, in particular for objects in the Galactic plane, they must generally be regarded as lower limits. A more reliable check of the X-ray determined $`N_\mathrm{H}`$ can come from the investigation of dust extinction maps. Using the afore mentioned Predehl & Schmitt (1995) relation, Verbunt at al. (1996) set an upper limit of $`2.2\times 10^{21}`$ cm-2 for the column density (i.e. in agreement with the DM value). From the more recent SFD map, as shown in Fig. 4, we note that PSR B1937+21 is almost at the center of an enhanced absorption region. The map has a pixel resolution of $`2\stackrel{}{.}37\times 2\stackrel{}{.}37`$ and spatial resolution of 6$`\stackrel{}{.}`$1 (FWHM). The enhanced region is then consistent with a point source at coordinates RA (J2000) $`19^\mathrm{h}`$ 39$`\stackrel{m}{.}`$5 Dec (J2000) $`21^{}`$ 37$`\stackrel{}{.}`$5. Given the mentioned spatial resolution, it is compatible with the position of PSR B1937+21. However at this position there is a strong radio and IR source: IRAS 19375+2130 ($``$ 4C21.53W $``$ G57$``$0.27 $``$ G57.55$``$0.27). This source was extensively investigated since the discovery of PSR B1937+21 (see e.g. Becker & Helfand 1983). It is listed as an ultra-compact H ii region and selected by Wood & Churchwell (1989) as a candidate site of OB star formation. To estimate the visual Galactic extinction toward the pulsar, we used the map value of neighboring regions: it is E(B$``$V) $`6`$ and a value of $`A_V19`$ is derived, giving $`N_\mathrm{H}=3.3\times 10^{22}`$ cm-2. This is some 50% higher than the X-ray derived $`N_\mathrm{H}`$ for PSR B1937+21. The observed high dust/gas ratio can then be justified without requiring a very large distance to the pulsar. However it has to be larger than the NE2001 model derived one (Cordes & Lazio 2002) as 3.6 kpc is a lower limit derived from timing (see below). Considering that the SFD map predictions are not very reliable at galactic latitudes $`|b|<5^{}`$ (Schlegel et al. 1998), the consistency with the X-ray derived value is remarkable.

### 5.2 The pulsar distance

Timing observations of the pulsar provide proper motion measurements and an upper limit on the timing parallax of $`\pi <0.28`$ mas (see Table 1). This information can be used to infer that the DM distance of 3.6 kpc is, in fact, a lower limit. An upper limit can be derived from the transversal velocity estimate given by the interstellar scintillation velocity $`v_{\mathrm{ISS}}50\pm 10`$ km s-1 found by Cordes (1990). Anyway, given the timing derived pulsar proper motion, the resulting limit of $`d=v/(4.74\mu )=v/(4.74\times 0.5)=21`$ kpc is not very constraining. While it has been shown that the scintillation velocity is indeed a good estimator for the true transversal velocity (Nicastro et al. 2001), for low velocities such as measured for MSPs (i.e. $`30`$ km s-1), systematic effects can result in large errors. A campaign of scintillation study could lead to a more reliable $`v_{\mathrm{ISS}}`$ estimate, but it is unlikely that a value lower than 30 km s-1 can be measured and therefore it can only be useful to set the Galactic upper limit distance of $`10`$ kpc. On the other hand, regular multi-frequency observations can help to remove the effects of “interstellar weather”, which degrades the timing precision of PSR B1937+21 significantly (e.g. Backer & Wong 1996). If these effects can be modeled, higher timing precision is obtainable in the future, possibly helping in setting tighter limits on the parallax measurements.

Heiles et al. (1983) studied the H i absorption spectra toward PSR B1937+21, 4C21.53W and another close-by source. By setting the Galactic tangent point to 5.4 kpc, a comparison of position and intensity of the lines for these three objects allowed them to conclude that 4C21.53W is at $`10.7`$ kpc and the pulsar at $`5`$ kpc from us. Frail & Weisberg (1990) revised this result using a new model for the rotation of the Galaxy. They adopt the tangent point distance as the pulsar distance lower limit: $`d>4.6\pm 1.9`$ kpc. Watson et al. (2003) have recently confirmed the position of the H ii region in the Perseus arm. Assuming a tangent point at a distance of 4.6 kpc, they derive a distance of 8.6 kpc by studying H2CO absorption and H110$`\alpha `$ emission lines (H110$`\alpha `$ is the 4874.1570 MHz radio recombination line). However the data are inconclusive on the position of the intervening molecular cloud, giving either 3.6 or 5.5 kpc as possible distances, for the two sides of the tangent point, respectively. With the pulsar at a distance $`>3.6`$ kpc, it is possible that the cloud contributes significantly to the absorption of the low energy X-rays.

We also investigated the 1.4 GHz NVSS data (Condon et al. 1998). Figure 5 shows the $`15^{}\times 15^{}`$ region around the pulsar position. PSR B1937+21 lies on the edge of a diffuse circular region (flux is $`50÷60`$ mJy/beam) and is moving away from it as indicated by the proper motion vector shown in Fig. 5. Incidentally, the vector orientation along the Galactic plane can explain why PSR B1937+21 has the lowest z-height of all MSPs. The strong source to the North is the aforementioned 4C21.53W and has a flux density of $`800`$ mJy, compared to the $`20`$ mJy of PSR B1937+21. The non-detection of the diffuse region in the VLA maps by Becker & Helfand (1983) can be explained by their use of a different array configuration. They used the B-array rather than the D-configuration as used for the NVSS, resulting in a much degraded sensitivity for diffuse emission. Though the circular shape of the region and the pulsar velocity vector could be suggestive of an association, the curved (externally “compressed”) shape of the H ii region indicates a relation with the latter (see Fig. 1 of Becker & Helfand 1983). A possible way to estimate the chance coincidence of PSR B1937+21 with the H ii “complex” is to select the sources of Watson et al. (2003) away from the inner arms, i.e. $`70^{}<l<50^{}`$ and $`|b|1^{}`$, and assume that they all have an associated “region” of $`4^{}`$ diameter (we recall that those H ii regions were selected from the IRAS *point source* catalogue, i.e. they have FWHM $`<2^{}`$ at 100 $`\mu `$m). The ratio of the total area of the selected 20 sources to $`20\times 2`$ square degrees is $`9\times 10^4`$. We note also that in the same area there are 15 known pulsars, two of which are MSPs444 http://www.atnf.csiro.au/research/pulsar/psrcat/. The coincidence at this galactic longitude of a peculiar pulsar with such a small region (at a distance of 8.6 kpc the region would be $`20`$ pc across, $`3.5`$ pc at 3.6 kpc) remains then remarkable.

### 5.3 Spectral index and magnetospheric emission

Table 3 lists the main emission properties of the three MSPs which show non-thermal emission only. In order to account for the contribution to the DC flux by the field source listed by Becker et al. (2003), but unresolved by the BeppoSAX MECS, we re-analysed the BeppoSAX observation of B1821$``$24. The pulsed fraction was compatible with 100%. This is consistent with the result by Danner et al. (1997) using ROSAT HRI data in the 0.1–2 keV band and explains the similar indices of the on-pulse BeppoSAX spectrum and the total CHANDRA spectrum, for which no pulse phase selection could be performed. As can be seen from Table 3, the “total” flux photon index of PSR B1937+21 is significantly softer than that measured for PSR J0218+4232 and PSR B1821$``$24. Also, we found that the spectrum of the pointlike DC component of PSR B1937+21 is softer than its pulsed spectrum. It is interesting to note that high-resolution CHANDRA observations of PSR J0218+4232 (Kuiper et al. 2002) confirmed the presence of a DC component of this source, earlier noted in Kuiper et al. (1998) and Mineo et al. (2000). The spectrum of the DC component of PSR J0218+4232 appeared to be softer than its pulsed component, like we found here for PSR B1937+21. Kuiper et al. (2000) concluded that we do not expect thermal emission from the surface of an old, recycled MSP due to cooling, and that the soft DC component is most likely due to reheating of the polar cap area by backflowing accelerated particles. This has to be investigated further for PSR B1937+21. A high spatial resolution observation may reveal field sources and/or diffuse emission also in the PSR B1937+21 field. Such information would allow us to better constrain both its total and pulsed emission properties as compared to those of PSR B1821$``$24 and PSR J0218+4232 for which higher resolution observations already exist. In spite of the large $`\dot{P}`$ difference between PSR B1937+21 and B1821$``$24, the X-ray production efficiencies of these two pulsars are similar; and similar are their $`B_L`$ and pulsed emission spectral indices. PSR J0218+4232 is a factor of 10 more efficient even though its $`B_L`$ is the weakest of the three; on the other hand it has the hardest spectrum. We finally note that using 4.6 kpc as the distance to PSR B1937+21 and as the distance to PSR J0218+4232 the 4 kpc upper limit recently set by Bassa et al. (2003), the “pulsed” X-ray luminosities of the three pulsars are almost identical. Identical are also the pulsed spectral indices if for PSR J0218+4232 we adopt the statistically more significant value of $`\alpha =1.14`$ from Kuiper et al. (2003).

All together these figures are not inconsistent with the claiming of a single $`L_X\dot{E}`$ law for MSPs and ordinary pulsars suggesting that a common magnetospheric emission process could be responsible for most, if not all, of the emitted X-ray flux in the band 2–10 keV. It is not clear if pulsar in 47 Tuc could be peculiar in this respect given that they have been studied only in the soft energy range 0.1–2.4 keV (Grindlay et al. 2002).

## 6 Conclusions

We detected the double peak profile of the fastest rotating pulsar known. The 1.3–10 keV pulsed fraction is 85% which becomes 100% in the 4–10 keV band. The DC component is significantly detected only in the low energy band 1.3–4 keV ($`46\pm 7\%`$ unpulsed photons). The pulse phase resolved spectral analysis confirms this result. The secondary (X-ray) peak is detected at high significance above 3–4 keV and the ratio primary/secondary decreases with energy. This suggests that the secondary peak has a harder spectrum. It was not possible to perform an absolute phase comparison with the radio profile, but phase separation comparison indicates, contrarily to the ASCA finding by Takahashi et al (2001), that the main radio and X-ray peaks are aligned. This indeed was confirmed by a Rossi-XTE observation which shows that the X-ray peaks are almost perfectly aligned with the radio giant pulses (Cusumano et al. 2003).

We measure a hard power-law spectral index $`\alpha 1.2`$ for the pulsed photons. This is similar to the value found for PSR B1821$``$24 and PSR J0218+4232 (Kuiper et al. 2003). We find a column density of $`N_\mathrm{H}2.1\times 10^{22}`$ cm-2 toward PSR B1937+21. It is marginally consistent with the integral Galactic value obtained by the Dickey & Lockman (1990) H i model of $`1.2\times 10^{22}`$ cm-2 but in agreement with the SFD map prediction, and gives a e-/$`N_\mathrm{H}`$ ratio of 1/100.

From the NVSS maps we note that the pulsar is located on the edge of a diffuse $`4^{}`$ circular region with the strong radio/IR source 4C21.53W and the pulsar located at opposite edges. The pulsar also seems to move away from this region with a transverse speed vector along the Galactic plane. Emission/absorption line observations have shown that 4C21.53W is very likely located in the Perseus arm at 8–9 kpc away from us (Watson et al. 2003) whereas it is unlikely the pulsar to be too far away from the Galactic tangent point, i.e. 4.6 kpc (Frail & Weisberg 1990). These results and the morphology of 4C21.53W/diffuse region suggest they are associated, with the pulsar aligned by chance.

A detailed study of the ISM toward the pulsar could possibly help modeling the (relatively) high timing noise of PSR B1937+21 and consequently further constrain the parallax mesurement. However, we consider $`5`$ kpc to be a good distance estimate.

A high spatial resolution CHANDRA X-ray observation can help both to verify the high $`N_\mathrm{H}`$ value and to study the continuum X-ray emission properties of its close surroundings. A better resolved 1.4 GHz radio map (obtainable from C configuration VLA archival data) can help to study the morphology of the “4C21.53 complex” and to confirm the chance alignment of the diffuse emission region with the pulsar.

###### Acknowledgements.

This research is supported by the Italian Space Agency (ASI) and Consiglio Nazionale delle Ricerche (CNR). BeppoSAX was a major program of ASI with participation of the Netherlands Agency for Aerospace Programs (NIVR). Its unburned pieces sinked in the Pacific Ocean on April 30, 2003, exactly on the date of its $`7^{\mathrm{th}}`$ birthday. We would like to thank the anonymous referee whose comments and suggestions improved the clarity of the paper.


format: text

text:
# Deriving the respiratory sinus arrhythmia from the heartbeat time series using Empirical Mode Decomposition

## 1 Introduction

Oscillations in heart rate on a beat by beat basis are a well-known phenomenon of great clinical relevance for both diagnostic and prognostic purposes. Many different factors (i.e. neural, chemical, hormonal modulations) contribute to the heart rate variability (HRV) . Among them, a major influence is exerted by the activity of the Sympathetic and Parasympathetic (vagal) branches of the Autonomic Nervous System (ANS) . In particular, respiration is a powerful modulator of HRV acting under the control of the ANS parasympathetic branch; in addition, abnormalities in respiratory modulation are often an indication of autonomic disfunction .

Traditionally HRV analysis is accomplished, in the frequency domain, via spectral analysis of the heartbeat time series (RR intervals) and, in the time domain, by extracting statistical indexes not related to specific cycle lengths. Parasympathetic activity is primarily reflected in the high-frequency (HF) component of the power spectrum (0.15-0.40 Hz) and it is related to the breathing frequency (respiratory sinus arrhythmia, RSA). More controversial is the interpretation of the low-frequency (LF) component (0.04-0.15 Hz), which is considered by some investigators as a marker of sympathetic modulation  and by others as influenced by both the sympathetic and parasympathetic activity . In the clinical environment a widely used index of sympathovagal balance is the ratio LF/HF.

For what concerns time domain indexes, they all measure high frequency variations in heart rate (and thus they are all highly correlated). An extensive introduction to HRV measurements can be found in .

Both time-domain and spectral methods share limitations imposed by the intrinsic nature of the RR series. In fact, changes in autonomic modulations can occur because of a variety of common factors such as changing position from sitting to standing, eating, running, sleeping, wakening. This can induce both a lack of stationarity in the series, compromising a correct application of these methods, and a shift of the LF and HF central frequencies, with a possible superposition of the two bands and a consequent misleading interpretation.

Other methods estimate RSA through a linear model by using an adaptive filter scheme having the respiratory signal as reference input . The first limitation of these approaches is the need of the respiratory signal. The second concerns the approximation of the linear model that is not always acceptable. Like almost all natural phenomena, HRV is the result of many nonlinearly interacting processes; therefore any linear analysis has the potential risk of underestimating, or even missing, a great amount of information content. On the other hand, complex nonlinear approaches  do not provide the tracking capability required in the high non-stationary context of RR interval series.

Recently the Empirical Mode Decomposition (EMD), a signal processing technique particularly suitable for nonlinear and nonstationary series, has been proposed  as a new tool for data analysis. This technique performs a time adaptive decomposition of a complex signal into elementary, almost orthogonal components that do not overlap in frequency. The extracted components have well behaved Hilbert transforms, from which the instantaneous frequencies can be calculated.

We applied EMD analysis to decompose the RR series into its components in order to identify, primarily, the respiratory oscillation. For comparison purpose, the respiratory signal was recorded simultaneously with the electrocardiogram (ECG). The existence of phase and frequency coupling between the RR component associated with breathing and the respiratory signal itself was then investigated.

## 2 Empirical Mode Decomposition

The empirical mode decomposition is a signal processing technique proposed in 1998 by Huang  to extract all the oscillatory modes embedded in a signal without any requirement of stationarity or linearity of the data. The goal of this procedure is to decompose a time series into components with well defined instantaneous frequency by empirically identifying the physical time scales intrinsic to the data, that is the time lapse between successive extrema.

Each characteristic oscillatory mode extracted, named Intrinsic Mode Function (IMF), satisfies the following properties: an IMF is symmetric, has a unique *local* frequency, and different IMFs do not exhibit the same frequency at the same time. In other words the IMFs are characterized by having the number of extrema and the number of zero crossings equal (or differing at most by one), and the mean value between the upper and lower envelope equal to zero at any point.

The algorithm operates through six steps:

1. Identification of all the extrema (maxima and minima) of the series $`X`$.
2. Generation of the upper and lower envelope via cubic spline interpolation among all the maxima and minima, respectively.
3. Point by point averaging of the two envelopes to compute a local mean series $`m`$.
4. Subtraction of $`m`$ from the data to obtain a IMF candidate $`h=Xm`$.
5. Check the properties of $`h`$:

* if $`h`$ is not a IMF (i.e. it does not satisfy the previously defined properties), replace $`X`$ with $`h`$ and repeat the procedure from step 1
    * if $`h`$ is a IMF, evaluate the residue $`r=Xh`$.
6. Repeat the procedure from step 1 to step 5 by *sifting* the residual signal. The sifting process ends when the residue $`r`$ satisfies a predefined stopping criterion.

At the end of the procedure we have a residue $`r`$ and a collection of $`n`$ IMFs, named $`h_i`$ ($`i=1,\mathrm{},n`$). The $`h_i`$ are generated being sorted in descending order of frequency and therefore $`h_1`$ is the one associated with the locally highest frequency. Furthermore the original $`X`$ can be exactly reconstructed by a linear superposition:

$$X=\underset{i=1}{\overset{n}{}}h_i+r.$$

Usually to eliminate *riding waves* and obtain symmetric waves, the sifting process has to be repeated a number of times. To verify whether $`h`$ owns the IMF properties (step 5 of the sifting procedure) we used a range-based criterion: $`h`$ is retained as IMF if the range of $`m`$ is a very low fraction of that of $`h`$ $`(<0.001)`$.

Analogously, for what concerns the stopping criterion of the EMD iterations, when all the interesting intrinsic modes of the original series have been extracted, we checked the range of $`r`$: the sifting process ends when the range of the residue is low with respect to that of the original signal ($`<10\%`$).

Figure 1 shows the starting point of a signal decomposition and the IMF candidate obtained after few iterations. Figure 2 shows an example of EMD decomposition of a simulated series obtained by linear composition of three different nonstationary oscillations.

The EMD procedure, according to the above specifications, was used for the decomposition of the RR interval series (tachograms) as described in the next section.

## 3 Experimental procedure

Thirteen healthy volunteers underwent an experimental session in which one electrocardiographic (ECG) derivation and the respiratory signal were simultaneously and continuously acquired at 1000 Hz sampling frequency. The ECG was recorded using standard electrodes, while the respiratory signal was detected through a polymeric piezoelectric dc-coupled transducer inserted into a belt wrapped around the chest. During the experimental session the subjects were comfortably sitting and breathing freely.

From the ECG signal, R waves, indicating the ventricular depolarization, were detected to obtain the tachogram, that is the series of the time intervals between the occurrence of two consecutive R peaks. The tachogram was then interpolated at 1 Hz, to obtain uniformly time spaced data.

The respiratory signals were bandpass filtered in the respiratory range 0.15-0.4 Hz  to remove artifacts mainly due to the positioning of the chest belt and its mechanical properties. The respiratory signal was then decimated according to the corresponding interpolated tachogram, thus making the two series syncronous. The tachograms of the 13 subjects were decomposed using EMD and each first component was analyzed together with the corresponding respiratory signal.

## 4 Instantaneous frequency analysis

The local frequency of each IMF can be expressed by the instantaneous frequency defined by virtue of the analytic signal of the original data.

Given a real function $`f(t)`$ with Fourier transform $`F(\nu )`$, the analytic signal $`f_a(t)`$ is defined as the inverse Fourier transform of the positive frequency part of $`F(\nu )`$:

$$f_a(t)=2_0^+\mathrm{}F(\nu )e^{j\mathrm{\hspace{0.17em}2}\pi \nu t}𝑑v.$$

The real part of $`f_a(t)`$ is just $`f(t)`$ while the imaginary part defines the Hilbert transform of $`f(t)`$ ($`f_H(t)`$). Thus $`f_a(t)=f(t)+jf_H(t)=a(t)\mathrm{exp}j\varphi (t)`$ in which:

$$a(t)=(f^2(t)+f_H^2(t))^{\frac{1}{2}}\varphi (t)=\mathrm{arctan}\frac{f_H(t)}{f(t)}.$$

The instantaneous frequency is defined as the rate of phase change:

$$\nu (t)=\frac{d\varphi (t)}{dt}.$$

We investigated the phase relationship between the first IMF and the respiratory signal by computing the instantaneous phases $`\varphi _1`$ and $`\varphi _2`$ from both signals and looking at their behavior by plotting $`\varphi _1`$ mod $`2\pi `$ versus $`\varphi _2`$, mod $`2\pi `$ (bi-plot).

When a phase locking exists, the points in the bi-plot are not uniformly distributed and some structure is visible . In case of linear phase coupling, the points in the bi-plot are ideally aligned along a straight line, in practice they cluster along a *stripe* whose slope indicates the frequency ratio (synchronization order) of the two signals.

The *stripe* slope was estimated through a generalized regression, that is detecting the direction of the eigenvector associated to the highest eigenvalue (maximum variance direction) via a Principal Components Analysis (PCA) procedure . The goodness of the relationship was measured by the residual variance, an index of the *spread* of points along the orthogonal direction. The instantaneous frequencies were computed using the toolbox *Time-frequency toolbox for Matlab* (TFTB) .

## 5 Results

The tachograms of all the 13 subjects were decomposed using EMD. The number of the oscillatory modes extracted ranged between 5 and 8 (median 7).

Since the breathing frequency is the highest frequency contributing to HRV, the first IMF (IMF1) generated was used for comparison with the recorded respiratory signal. For each subject IMF1 and the recorded respiratory signal were Hilbert-transformed and the instantaneous phases computed. Figure 3 shows the tachogram and and the simultaneous respiratory signal of one subject.

The relative phase distribution $`\mathrm{\Delta }\varphi `$ = $`\varphi _1\varphi _2`$ shown in figure 4 exhibits a sharp peak suggestive of a coupling between $`\varphi _1`$ and $`\varphi _2`$. An example of the close correspondence between the two series of instantaneous frequencies is presented in figure 5.

The average frequency ratio between each couple of signals is $`1.015\pm 0.017`$ which is very close to 1, as expected for a synchronization of order 1:1. The average residual variance is 3.7%, indicating a good linear relationship.

## 6 Discussion

The Empirical Mode Decomposition technique allows the analysis of complex, nonlinear and nonstationary time series, through their decomposition into a limited number of *elementary* modes having interesting local properties. In the EMD applications reported in literature , the extracted modes are speculatively associated with specific physical or physiological aspects of the phenomenon investigated. In this application we demonstrated the association of the first IMF extracted from a tachogram with the simultaneously recorded respiratory signal.

The EMD implementation required some computational skills essentially to cope with the selection of the stopping criteria for both sifting and EMD processes and the management of the end points for cubic splines interpolation. The selected criteria were not necessarily unique: in any case the procedure is robust with respect to a variety of different solutions.

The peculiar property of each IMF of having a single local frequency was particularly suitable to compute the instantaneous phases and frequencies after Hilbert transformation. According to , in order to verify the synchronization between the two series, we focused on the phases relationship. We preferred to work with instantaneous phases instead of frequencies for that the former resulted to be more stable in time. The series of frequencies shown in figure 6 are indeed not instantaneously derived, but computed with a smoothing factor from the toolbox TFTB.

The linear relationship and the generalized regression slope closest to 1 confirmed the synchronization between $`\varphi _1`$ and $`\varphi _2`$ and therefore the substantial correspondence between the IMF1 and respiratory time series. The importance of this finding is essentially two-fold for the physiologist. First, the instantaneous respiratory contribution can now be available even in all those retrospective studies where the breathing signal was lacking. Second, the possibility of selectively and efficiently eliminate the respiratory component from the tachogram makes possible a more accurate evaluation of the LF role in ANS studies.

Our goal is now to proceed with this approach in the attempt of attributing a physiological meaning to the other modes embedded in the tachogram in order to improve the knowledge of the pathophysiological mechanisms underlying heart rate variability.


text:
# Hadronic transitions from the lattice

## 1 INTRODUCTION

Rather few hadronic states are actually stable under strong interactions. Because of this, a full understanding of hadrons from QCD will involve developing techniques to treat unstable hadrons and to evaluate their hadronic transition strengths.

An area of hadronic physics that is of considerable interest is that of gluonic excitations: the exploration of hadrons with a non-trivial gluonic content. The most clear case is that of glueballs, however, any glueball states will mix with quark-antiquark mesons in practice. This mixing can be thought of as mediated by a hadronic transition. Moreover, the eventual hadronic decays of such mixed mesons also need to be evaluated to give a clearer picture of the relationship between the physical spectrum and the ingredients (glueballs and quark-antiquark mesons).

The other natural area to explore gluonic excitations is that of hybrid mesons. These are quark-antiquark mesons with an additional non-trivial gluonic component. The benchmark example is a spin-exotic meson (so-called because it cannot be made from combining a quark and antiquark). These spin-exotic hybrid mesons will also be unstable, and study of their decay transitions will guide the interpretation of experimental candidates.

An understanding of methods to treat purely hadronic transitions will also help in developing methods to explore weak decays to hadronic final states (e.g. $`K\pi \pi `$).

The lattice QCD approach is the only quantitative non-perturbative method but the formulation in euclidian time causes problems with the treatment of decays as I now discuss.

### 1.1 Euclidian time

Hadron masses are extracted from lattice QCD calculations using two point correlators, however, their behaviour as $`e^{Mt}`$ in euclidian time will not be appropriate for hadrons that can decay. The most naive modification of the lattice QCD formalism caused by the introduction of decay widths is the replacement $`MMi\mathrm{\Gamma }/2`$ where $`\mathrm{\Gamma }`$ is the decay width. This gives a behaviour as $`e^{Mt}e^{i\mathrm{\Gamma }t/2}`$ but this is not consistent: the correlation must be positive definite yet this expression oscillates. Here I discuss this apparent paradox more fully .

The euclidian time propagator of an unstable particle which can decay into two stable particles of mass $`\mu `$ can be expressed as

$$G(t)=\frac{1}{\pi }_{2\mu }^{\mathrm{}}𝑑Ee^{Et}\rho (E)$$
(1)
where we can express the spectral density for a resonance of mass $`M`$ as

$$\rho (E)=\frac{1}{2M}\frac{\mathrm{\Gamma }(E)/2}{(ME)^2+(\mathrm{\Gamma }(E)/2)^2}$$
(2)
As an illustration, consider taking $`\mathrm{\Gamma }(E)`$ as a constant with $`\mathrm{\Gamma }(M2\mu )`$ and $`(M2\mu )t1`$. Then we have, in this approximation,

$$2MG(t)=e^{Mt}\mathrm{cos}(\mathrm{\Gamma }t/2)+\frac{e^{2\mu t}\mathrm{\Gamma }}{2\pi (M2\mu )^2t}$$
(3)
This expression, when evaluated, shows that the threshold contribution (the second term) dominates at larger $`t`$ and hence the oscillating behaviour of the first term will be obliterated. This clarifies the apparent paradox mentioned above concerning the sign of $`G(t)`$: the contribution from the lower lying (in energy) component of the resonance is enhanced in euclidian time. Indeed this can be seen directly from eq. 1.

By contrast, the Minkowski behaviour of the two particle correlator under the same assumptions is

$$2MG_M(t)=e^{iMt}e^{\mathrm{\Gamma }t/2}\frac{ie^{2i\mu t}\mathrm{\Gamma }}{2\pi (M2\mu )^2t}$$
(4)
for which the second term is almost always negligible.

In principle, a lattice determination of the spectral function $`\rho (E)`$ from measurements of $`G(t)`$ using eq. 1 would be sufficient to give a full description of hadronic decays. This is not a simple task: the inverse Laplace transform need to obtain $`\rho (E)`$ from $`G(t)`$ is not numerically stable. It has to be stabilised by making model assumptions about $`\rho (E)`$, for example as discussed in refs. . Another possibility is that the maximum entropy method may provide sufficient numerical stability and this has been tested in model cases .

The same observation that the region near to threshold dominates is at the root of the conclusions  of Maiani and Testa that two body states have unappetizing properties in euclidian time.

### 1.2 Particles in a box

In practice this issue is less relevant since the lattice evaluations are performed in a finite spatial volume with periodic boundary conditions. This finite spatial size of the lattice implies that two-body states are actually discrete. By measuring their energy very precisely as the spatial volume is varied, it is possible  to extract the scattering phase shifts and hence decay properties.

Consider the two body states of two pions of mass $`\mu `$ with momentum $`\pm 𝐤`$ which will have total energy $`E=2(\mu ^2+𝐤^2)^{1/2}`$. Then if they are non-interacting, their momentum on a spatial hyper-torus of size $`L`$ will be discrete with $`𝐤=2\pi 𝐧/L`$ where $`𝐧`$ has integer components (and for future use we define $`N`$ such that there are $`N`$ values of $`𝐧`$ having a given $`n=|𝐧|`$, so $`N=6`$ if $`n=1`$). If the two pions interact with a finite range of interaction, compared to the spatial volume, this implies that the finite size effect will result in an momentum shift of order $`L^3`$, provided that $`L`$ is large enough to contain one pion without undue distortion. The detailed analysis of Lüscher  expresses the shift in this momentum value when there are interactions among the pions with phase-shift $`\delta (k)`$, provided $`E`$ is below the inelastic threshold, as

$$d(E^2)=4d(k^2)=\frac{16\pi N\mathrm{tan}\delta (k)}{L^3k}$$
(5)
to leading order in $`L^1`$, and indeed a much more precise expression with sub-leading terms is available through Lüscher’s work .

Now when there is a resonance at $`m_R`$ near to a two particle state, there will be an opportunity to deduce the resonance width. Then, one can describe the phase shift $`\delta (k)`$ by using the expression for elastic $`\pi \pi `$ scattering dominated by this resonance pole:

$$\mathrm{tan}\delta (k)=\frac{\mathrm{\Gamma }(k)}{2(m_RE)}$$
(6)
so that determinations of $`\delta (k)`$ at two or more values of $`k`$ will enable the resonance parameters $`m_R`$ and $`\mathrm{\Gamma }(k_R)`$ to be evaluated, where $`m_R^2=4(\mu ^2+k_R^2)`$. The most promising way to determine the two particle energy accurately is by using several lattice operators (including some built like two pions and some built like the resonance under study) as sources and sinks to get the best determination of the energy levels from fitting the matrix of correlators. Note that in practice it will be difficult to obtain sufficiently accurate energy determinations of the lowest two-particle state and it will be even more difficult to determine the energy level of the next heavier state with larger momenta.

### 1.3 A direct approach

It should be possible to extract the hadronic transition amplitude directly from the lattice (rather than via energy level shifts), and this we now discuss.

Using suitable lattice operators to create $`R`$ at $`t=0`$ and annihilate a two-pion state with momenta $`k`$ and $`k`$ at time $`t`$, the contribution to the correlator from a $`R`$ state with mass $`m_R`$ and a two-body state with energy $`E_{\pi \pi }`$ is given to leading order in $`x`$ by

$$C_{R\pi \pi }(t)=\underset{t_1}{}he^{m_Rt_1}xe^{E_{\pi \pi }(tt_1)}b$$
(7)
where there is a summation over the intermediate $`t`$-value $`t_1`$ and where $`h`$ and $`b`$ are the amplitudes to make each state from the lattice operators used and $`x`$ is the required transition amplitude $`R|\pi \pi `$. Here we are assuming that the states $`R`$ and $`\pi \pi `$ are normalised to unity. By obtaining $`h`$ and $`b`$ from the $`RR`$ and $`\pi \pi \pi \pi `$ correlators, one can attempt to extract $`x`$.

The complication, however, is that removal of excited state contributions is tricky, even if variational methods are used to construct improved lattice operators to create the ground states. For example, if $`m_RE_{\pi \pi }>0`$ then the transition time $`t_1`$ will be preferentially near 0 (since the heavier state then propagates less far in time) and one can complete the sum over $`t_1`$ obtaining a $`t`$-dependence of eq. 7 as $`e^{E_{\pi \pi }t}`$. This same $`t`$-dependence would be obtained if the state with mass $`m_R`$ were to be replaced with an excited state with an even heavier mass. Thus one cannot separate the ground state and excited state contributions even in principle. See ref.  for a fuller discussion.

The way forward is that if $`m_R=E_{\pi \pi }`$, the ground state contributions to eq. 7 have a $`t`$-dependence as $`te^{E_{\pi \pi }t}`$ whereas any excited state contributions behave as $`e^{E_{\pi \pi }t}`$ as above. So now we do have a way to isolate the required ground state contribution:

$$xt=\frac{C_{R\pi \pi }(t)}{[C_{RR}(t)C_{\pi \pi \pi \pi }(t)]^{1/2}}+𝒪(x)+𝒪(x^3t^3)$$
(8)
Note that this separation is only by a power of $`t`$ which is less than the case for two-point correlations where the excited state contributions are suppressed by an exponential $`e^{(m_R^{}m_R)t}`$.

In practice the requirement of energy equality can be relaxed. Defining $`\mathrm{\Delta }=m_RE_{\pi \pi }`$, then the ground state contribution to the expression of eq. 8 evaluates to $`2x\mathrm{sinh}(t\mathrm{\Delta }/2)/(t\mathrm{\Delta })=x(1+(t\mathrm{\Delta })^2/24+\mathrm{})`$. So this will be equivalent to the expression with $`\mathrm{\Delta }=0`$ provided

$$(m_RE_{\pi \pi })t5$$
(9)
So far we have described the behaviour of the $`C_{R\pi \pi }(t)`$ in the limit of small $`x`$. The first correction term from multiple transitions (to eq. 8) will be of order $`x^3t^3/6`$ so we need $`xt1`$. As well as these corrections which are of higher order in $`xt`$, one must also consider the intrinsic mixing of the initial $`R`$ state with $`\pi \pi `$ (and vice versa). This intrinsic mixing (i.e. coming from the lattice operators used to create the states having an admixture - see fig. 1 for an example) is expected be of order $`x/`$ where $``$ is the energy of the quark pair and so will contribute a term like $`xe^{Et}/`$. This is a contribution similar to that from excited states and so will be dominated at large $`t`$ by the $`xte^{Et}`$ term we are looking for. So we need both $`xt`$ to be small and $`t`$ to be large. This implies that $`x`$ must be small for this direct approach. And indeed, it is only when $`x`$ is small that the more rigorous approach of determining the energy shifts will be computationally difficult to implement.

In summary, when $`t(m_RE_{\pi \pi })5`$; $`xt1`$ and $`(E^{}E)t1`$ (with $`E^{}E`$ the energy gap to the first excited state), it is possible to estimate the hadronic transition amplitude directly, and this we will explore in detail for $`\rho `$ decay.

### 1.4 Comparison

In order to show the relationship of this direct approach with Lüscher’s approach, consider the case when a two particle state is close in energy to the resonance. It should be possible to arrange that energy levels are sufficiently close by choosing the lattice spatial length $`L`$. Then once $`x`$ has been determined in this situation where the resonance does not decay significantly (because the two particle level has the same energy as the resonance), one can assume that the coupling strength determined can be used in a larger spatial volume where decay is energetically allowed.

The most direct way to evaluate resonance decay is using first order perturbation theory (Fermi’s Golden Rule) which implies a transition rate

$$\mathrm{\Gamma }=2\pi x^2\rho (E)$$
(10)
where the angle brackets indicate that an average over spatial directions will be needed. For a decay from the centre of mass with relative momentum $`k`$, the density of states $`\rho (E)=L^3kE/(8\pi ^2)`$.

Now to compare with Lüscher’s approach, we need to evaluate the energy shift of the two particle state in a finite box (which is exactly $`N`$fold degenerate). This can be estimated from the mixing of the two energy levels which are close to each other, using second order perturbation theory (in $`x`$) which gives an energy shift, on the lattice, of

$$dE=\frac{Nx^2}{m_RE},$$
(11)
Thus when the resonance lies above the two particle state, the mixing will move the resonance mass up and the two particle energy down, in the usual way as mixing causes repulsion of levels.

Using the relationship between $`x`$ and $`\mathrm{\Gamma }`$ from eq. 10, we obtain

$$dE=\frac{4\pi N\mathrm{\Gamma }}{(m_RE)kEL^3}$$
(12)
which is exactly the same expression as would be obtained at leading order in $`L^1`$ from Lüscher’s formalism as presented above. This is as it should be, but does illustrate that an estimate of $`x`$ directly from the lattice may provide a useful way to evaluate these decay effects.

Indeed the mass shifts are of second order in $`x`$ so will be small in general, and hence difficult to determine with sufficient precision. The direct measurement of the transition giving $`x`$ from the lattice is thus an attractive prospect.

## 2 VECTOR MESON DECAY

### 2.1 Transitions

For on-shell transitions, it is possible to estimate hadronic transition strengths as described above and here we explore this approach further for the case of $`\rho `$ meson decay to $`\pi \pi `$, following ref. .

The situation we shall analyse is represented by the energy spectrum shown in fig. 2, here neglecting any interactions among the states. Note that no decay can actually take place with these lattice parameters. This is , however, optimum for our analysis since when the energy levels are approximately degenerate, the hadronic transition amplitude can be evaluated more effectively.

We evaluate correlations between lattice operators creating both a $`\rho `$ meson and a $`\pi \pi `$ state, using a stochastic method with sources on time-slices  to evaluate the quark diagrams shown in fig. 3 from 20 gauge configurations  with $`N_f=2`$ sea quarks of mass corresponding to about 2/3 of the strange-quark mass. These correlations, normalised by the two point functions as appropriate, are illustrated in fig. 4. Here the off-diagonal case (labelled $`\rho _1\pi _1\pi _0`$) shows the important feature that it grows approximately linearly with increasing $`t`$. As emphasized above, this linear growth will only occur for on-shell transitions and this is essentially the case here.

To extract an estimate of the transition amplitude, we wish to evaluate $`x=\rho |\pi \pi `$ where these states are normalised on the lattice (to unity). If higher excited states are neglected, one can make a two state model (i.e. $`\rho `$ and lightest $`\pi \pi `$ state) with this transition amplitude and evaluate the contributions as shown by the curves in fig. 4. Because the location (at $`t_1`$ between 0 and $`t`$) of any transition is not known, one has potentially very serious contributions from excited states. As discussed above, these excited state contributions can be avoided if the transition is approximately on-shell, when the linear dependence on $`t`$ of the off-diagonal transition (with value $`xt`$) is a unique signature of this on-shell transition that we wish to extract.

Note that for the case of $`\rho _0\pi _1\pi _1`$, the on-shell condition is much less well satisfied but the relative momentum of the pions in the centre of mass is twice as large as for $`\rho _1\pi _0\pi _1`$ and hence $`x`$ should be approximately twice as large (as we find), since for a P-wave decay there will be a momentum factor in the transition amplitude.

A further check of the extraction of $`x`$ comes from the box diagram, fig. 3d, which will have a contribution behaving as $`x^2t^2/2`$ arising from a $`\rho `$ intermediate state - see the data labelled $`\pi _1\pi _0\pi _1\pi _0`$ in fig. 4. Again this contribution does seem to be present at the required level.

Thus we have qualitative agreement from the three-point evaluations with two different relative momenta and from the four-point evaluation that $`xa0.06`$.

A more quantitative estimate of $`x`$ can be made  by reducing excited state contributions as much as possible. Thus one uses fits to remove excited state contributions from the normalising two-point functions and then evaluates the slope of the ratio of eq. 8 versus $`t`$ to remove the $`𝒪(x)`$ term also.

Note that the above methods can be used in the quenched approximation  to estimate hadronic transition amplitudes. They depend on assuming that $`xa`$ is fairly weak, as we indeed find.

### 2.2 Energy shifts

A more rigorous approach is to focus on energy values. When two levels are close (our on-shell condition), then they will mix and the resultant energy shifts give relevant information  as discussed above. Moreover we can estimate these energy shifts from our $`x`$-value which provides more cross checks. These energy shifts can only be studied using dynamical fermions.

From a full analysis (i.e. measuring all the quark diagrams illustrated in fig. 3 and fitting the matrix of correlators) we obtain the energy shift of the $`\pi _1\pi _0`$ state (i.e. the un-binding energy), as needed in Lüscher’s approach, as 0.02(2) (in lattice units) upward which is consistent but not sufficiently accurate to be of significant use. The energy shift of the $`\rho _1`$ state can, however, be determined in this case because of a lattice artifact. The $`\rho `$ with momentum 1 (in lattice units of $`2\pi /L`$) can have its spin aligned parallel to the momentum axis (P) or perpendicular to it (A). Because the $`\pi _0\pi _1`$ state has relative momentum along a lattice axis and the transition from $`\rho `$ to $`\pi \pi `$ has orbital angular momentum L=1 (so a distribution like $`\mathrm{cos}\theta `$ in the centre of mass), only the parallel state (P) can mix with this two pion state. This situation is not typical - it occurs because the $`\rho `$ with momentum 1 only has a transition to a $`\pi _0\pi _1`$ state with a pion having momentum in the same direction. This mixing will not be present in the quenched approximation, so this provides a direct opportunity to see the effect of the two pion channel on the $`\rho `$ in unquenched studies. We do indeed find  such a mass splitting between the P and A orientations of the $`\rho `$ due to mixing for $`N_f=2`$ but not for $`N_f=0`$ as shown in fig. 5. Moreover the magnitude of this energy shift ($`0.026\pm 7`$ in lattice units) is consistent with other determinations of the transition strength.

### 2.3 Phenomenology

The basic assumption is that the transition from $`\rho `$ to $`\pi \pi `$ is given by an effective interaction with a finite spatial extent, this is usually summarised by an effective lagrangian where we normalise the coupling as $`\overline{g}^2=\mathrm{\Gamma }ME/k^3`$ in terms of the decay width. Then, provided the lattice spatial size is big enough that the hadrons are not distorted, our lattice situation (where no decay occurs) can be used to determine $`\overline{g}`$ and this can then be used to predict decay widths when the lattice volume is increased so that the minimum momentum becomes small. In our case we will also need to reduce the quark mass to allow decay since we have $`m_\pi /m_\rho =0.58`$.

One complication is that the transition which is closest to being on-shell is $`\rho _1\pi _1\pi _0`$ which is not in the centre of mass. A generalisation of Lüscher’s approach has been made  which allows for this. The more direct method we have described here can also be extended to this case . There are possible problems in this treatment of decays of moving particles. For instance, on the lattice a boost to bring the $`\rho `$ to rest will not bring the two pion state to have zero net momentum also. This is because the lattice situation allows energy non-conservation which, upon boosting, will imply some momentum non-conservation.

From our lattice studies, we deduce that $`ax=0.06_1^{+2}`$ for $`\rho _1\pi _1\pi _0`$. Translating  this lattice transition amplitude to the continuum normalisation, gives $`\overline{g}=1.40_{23}^{+27}`$. Using the observed $`\rho _1`$ energy shift gives another estimate, namely $`\overline{g}=1.56_{13}^{+21}`$. These two values agree well, which is good support for our scheme of extracting hadronic transition strengths.

Note that our lattice values would need to be extrapolated to light sea-quarks and to the continuum limit to allow all sources of systematic error to be explored in the comparison with the experimental situation. Nevertheless, these two values agree well with the values extracted from experimental data on decays of $`\rho `$, $`K^{}`$ and $`\varphi `$ mesons, namely $`\overline{g}1.5`$. This further underlies the viability of our methods.

## 3 HYBRID MESON DECAY

Hybrid mesons are those with non-trivial gluonic excitations. The spectrum of such mesons has been determined from lattice studies, both for the case of heavy quarks and light quarks. Hybrid mesons have allowed couplings to meson meson systems. If these coupling are large, then the hybrid mesons decay widths will be big and such mesons will be difficult to identify experimentally. Here we present recent lattice results on these hadronic transitions.

### 3.1 Hybrid mesons with heavy quarks

One well defined way to treat heavy quarks (such as the $`b`$-quark) on a lattice is to use the leading order of the HQET which corresponds to the static approximation. Here we take the heavy quarks to be fixed at locations $`R`$ apart in the $`z`$-direction.

Consider $`Q\overline{Q}`$ states with static quarks in which the gluonic contribution may be excited. We classify the gluonic fields according to the symmetries of the system. This discussion is very similar to the description of electron wave functions in diatomic molecules. The symmetries are (i) rotation around the separation axis $`z`$ with representations labelled by $`J_z`$ (ii) CP with representations labelled by $`g(+)`$ and $`u()`$ and (iii) C$`\mathrm{}`$. Here C interchanges $`Q`$ and $`\overline{Q}`$, P is parity and $`\mathrm{}`$ is a rotation of $`180^0`$ about the mid-point around the $`y`$ axis. The C$`\mathrm{}`$ operation is only relevant to classify states with $`J_z=0`$. The convention is to label states of $`J_z=0,1,2`$ by $`\mathrm{\Sigma },\mathrm{\Pi },\mathrm{\Delta }`$ respectively. The ground state ($`\mathrm{\Sigma }_g^+`$) will have $`J_z=0`$ and $`CP=+`$.

The exploration of the energy levels of other representations has a long history in lattice studies . The first excited state is found to be the $`\mathrm{\Pi }_u`$. This can be visualised as the symmetry of a string bowed out in the $`x`$ direction minus the same deflection in the $`x`$ direction (plus another component of the two-dimensional representation with the transverse direction $`x`$ replaced by $`y`$), corresponding to flux states from a lattice operator which is the difference of U-shaped paths from quark to antiquark of the form $``$.

Recent lattice studies  have used an asymmetric space/time spacing which enables excited states to be determined comprehensively. These results confirm the finding that the $`\mathrm{\Pi }_u`$ excitation is the lowest lying and hence of most relevance to spectroscopy.

From the potential corresponding to these excited gluonic states, one can determine the spectrum of hybrid quarkonia using the Schrödinger equation in the Born-Oppenheimer approximation. This approximation will be good if the heavy quarks move very little in the time it takes for the potential between them to become established. More quantitatively, we require that the potential energy of gluonic excitation is much larger than the typical energy of orbital or radial excitation. This is indeed the case , especially for $`b`$ quarks. Another nice feature of this approach is that the self energy of the static sources cancels in the energy difference between this hybrid state and the $`Q\overline{Q}`$ states. Thus the lattice approach gives directly the excitation energy of each gluonic excitation.

The $`\mathrm{\Pi }_u`$ symmetry state corresponds to excitations of the gluonic field in quarkonium called magnetic (with $`L^{PC}=1^+`$) and pseudo-electric (with $`1^+`$) in contrast to the usual P-wave orbital excitation which has $`L^{PC}=1^{}`$. Thus we expect different quantum number assignments from those of the gluonic ground state. Indeed combining with the heavy quark spins, we get a degenerate set of 8 states:

Note that of these, $`J^{PC}=1^+,0^+`$ and $`2^+`$ are spin-exotic and hence will not mix with $`Q\overline{Q}`$ states. They thus form a very attractive goal for experimental searches for hybrid mesons.

### 3.2 Hybrid decays

Within this static quark framework, one can explore the decay mechanisms. One special feature is that the symmetries of the quark and colour fields about the static quarks must be preserved exactly in decay, hence the light quark-antiquark pair produced must respect these symmetries. This has the consequence that the decay from a $`\mathrm{\Pi }_u`$ hybrid state to the open-$`b`$ mesons ($`B\overline{B},B^{}\overline{B},B\overline{B^{}},B^{}\overline{B^{}}`$) will be forbidden  if the light quarks in the $`B`$ and $`B^{}`$ mesons are in an S-wave relative to the heavy quark (since the final state will have the light quarks in either a triplet with the wrong $`CP`$ or a singlet with the wrong $`J_z`$). The decay to $`B^{}`$-mesons with light quarks in a P-wave is allowed by symmetry but not energetically.

In the heavy quark limit, the only allowed decays are when the hybrid state de-excites to a non-hybrid state with the emission of a light quark-antiquark pair. Since the $`\mathrm{\Pi }_u`$ hybrid state has the heavy quark-antiquark in a triplet P-wave state, the resulting non-hybrid state must also be in a triplet P-wave since the heavy quarks do not change their state in the limit of very heavy quarks. Thus the decay for $`b`$ quarks will be to $`\chi _b+M`$ where $`M`$ is a light quark-antiquark meson in a flavour singlet. This proceeds by a disconnected light quark diagram and it would be expected  that the scalar or pseudoscalar meson channels are the most important (i.e. they have the largest relative OZI-rule violating contributions). This transition can be estimated on a lattice when the initial and final energies are similar. This is the case  for the $`\mathrm{\Pi }_u`$ de-excitation to ground state gluonic field plus $`f_0`$ meson when the interquark separation is around 0.2 fm which allows a lattice evaluation of the hadronic transition strength - see fig. 6.

Note that in fig. 6 for $`R=0.1`$ fm the $`\mathrm{\Pi }_u`$ state is unstable to decay and, since we are using $`N_f=2`$ flavours of sea-quark, this decay is enabled and we might expect to see the lower energy (i.e. $`\mathrm{\Sigma }_g^++f_0`$) instead for this energy level. In practice, as we shall find later, the transition strength is quite small with $`x0.01`$ so the admixture of the lighter state is very small indeed (it is of order $`x^2`$) and hence it will only be significant at very large $`t`$-values, much larger than those used here.

Indeed the dominant transition is found to be with $`M`$ as a scalar meson and the signal for the hadronic transition strength $`x`$ is shown in fig. 7. Making use of the transition strength evaluated for a range of $`R`$-values around 0.2 fm and folding in the wavefunction overlaps, we can estimate the decay width of the hybrid spin-exotic meson $`H`$ (e.g. $`J^{PC}=1^+`$) which has $`\mathrm{\Pi }_u`$ excited glue. This corresponds to a width of around 100 MeV for the transition $`H\chi _b+f_0`$, whereas when $`M`$ is an $`\eta `$ or $`\eta ^{}`$ meson the transition strength is found to be less than a few MeV. There will be modifications to this analysis coming from corrections to the heavy quark limit (of order $`1/m_Q`$ where $`m_Q`$ is the heavy quark mass) which might allow hybrid meson transitions to $`B\overline{B}`$, etc, but these have not been evaluated yet.

In this heavy quark (or static) limit, the spin-exotic and non spin-exotic hybrid mesons are degenerate. For the latter, however, the interpretation of any observed states is less clear cut, since they could be conventional quark antiquark states. Moreover, the non spin-exotic hybrid mesons can mix directly (i.e. without emission of any meson $`M`$) with conventional quark antiquark states once one takes into account corrections (of order $`1/M_Q`$) to the static approximation applicable for heavy quarks with physical masses.

It is encouraging that the decay width comes out as relatively small, so that the spin-exotic hybrid states should show up experimentally as sufficiently narrow resonances to be detectable. This decay analysis does not take into account heavy quark motion or spin-flip and these effects will be significantly more important for charm quarks than for $`b`$-quarks.

## 4 CONCLUSIONS

It is possible to extract information about hadronic decays from euclidian lattice studies. The most rigorous method is to measure the small energy shifts in the two body states precisely as the lattice volume is varied. This then allows the phase shift to be extracted which gives information on resonance properties. This energy shift is only enabled in dynamical fermion simulations if quark-antiquark production is needed in the decay. This makes this a very difficult route to follow in practice.

We have also sketched an alternative which is less rigorous but will enable estimates to be obtained more readily. This direct method was tested in the case of $`\rho \pi \pi `$ where it gave consistent results and, moreover, results in agreement with experiment. This is a good indication that the method is reliable in practice. The method works best when the lattice transition amplitude ,$`x`$, is relatively small and it is encouraging that the method gives good indications when $`ax0.06`$. We then presented an application to hybrid meson decays: a case where the experimental result is not known, but where the magnitude of $`x`$ is indeed small. Our result is of experimental relevance: the knowledge of the magnitude of the decay width confirms that these hybrid states are likely to be accessible to straightforward study.

## 5 ACKNOWLEDGEMENT

I have benefited from helpful discussions with my colleagues in the UKQCD collaboration, especially Craig McNeile.