text:
# Multisite observations of the PMS 𝛿 Scuti star V351 OriBased on observations collected at the: Loiano Observatory (Italy), San Pedro Martir Observatory (México), OGS at Teide Observatory (Spain), JKT at Roque de los muchachos Observatory (Spain), Beijing Astronomical Observatory (China), SARA Observatory at Kitt Peak, which is owned and operated by the Southeastern Association for Research in Astronomy.

## 1 Introduction

$`\delta `$ Scuti stars are usually associated with the main-sequence (MS) or post-MS evolutionary phases. To date, about ten intermediate-mass pre–main-sequence (PMS) stars have been found to pulsate with timescales typical of $`\delta `$ Scuti variables and a growing interest in these young objects has developed in the last few years (see Marconi et al. 2001, hereafter M01; Marconi, Palla & Ripepi 2002, hereafter M02, and references therein). The first two candidates were discovered by Breger (1972) in the young cluster NGC 2264. The existence of PMS $`\delta `$ Scuti stars was later confirmed by observations of the short term variability of the Herbig Ae stars HR 5999 (Kurtz & Marang 1995) and HD 104237 (Donati et al. 1997).

The initial evidence has stimulated the theoretical investigation of the PMS instability strip, based on nonlinear convective hydrodynamical models (Marconi & Palla 1998). As a result, the topology of the instability strip for the first three radial modes was identified and a list of possible PMS pulsating candidates, mainly Herbig Ae stars, was selected on the basis of their spectral types. Then, a series of observing runs has been carried out searching for $`\delta `$ Scuti-type photometric variations with periods of minutes to a few hours and amplitudes less than a few tenths of a magnitude among those stars inside or close to the theoretical instability strip. These studies have confirmed the $`\delta `$ Scuti variability of HD 104237 (Kurtz & Muller 1999), HR 5999 (Kurtz & Catala 2001), and six additional candidates have been identified (see Kurtz & Muller 2001; M01; M02).

Among the new discoveries, V351 Ori is a particularly interesting object. The star is characterized by strong long-term photometric variations, attributed to extinction by circumstellar dust clouds. This property has been considered a sign of youth, signaling the transition from an active PMS phase to the arrival on the main sequence. Similar behavior has been observed in other well known Herbig Ae stars, such as BN Ori, LkH$`\alpha `$ 234 and UX Ori (e.g., Herbst & Shevchenko 1998). The spectral properties of V351 Ori include the presence of H$`\alpha `$ emission and a pronounced inverse P Cygni profile (van den Ancker et al. 1996), the SI $`\lambda `$1296 Å and OI $`\lambda `$1304 Å ultraviolet lines (Valenti et al. 2000), and a broad infrared excess (van den Ancker et al. 1997). Altogether, these features clearly reveal an active phase of interaction with the circumstellar environment. On the other hand, the possibility that V351 Ori is a post-MS star has been put forward in the literature. Evidence in this sense comes from the lack of a reflection nebulosity, the absence of dense molecular gas around the star, and a space velocity which indicates that the star is approaching Orion and hence is not kinematically linked to the star forming region. Spectral analysis of hydrogen absorption lines also show some peculiarities (Koval’chuck & Pugach 1998; see also M01 for details). Recently, on the basis of high dispersion echelle spectra, Balona, Koen & van Wyk (2002, hereafter BKW) have provided further evidence that V351 Ori is indeed a PMS star.

As already remarked, the interest for V351 Ori lies in the fact that it pulsates. Marconi et al. (2000) and M01 have interpreted the frequencies observed in the light curve of V351 Ori as due to pulsations in a mixture of several radial modes. In particular, in M01 four of the six observed frequencies have been simultaneously reproduced by a linear nonadiabatic pulsation model with mass $`M_{}=1.8M_{}`$, luminosity $`L_{}=13.87L_{}`$ and effective temperature $`T_{\mathrm{eff}}=7350K`$. However, the comparison of the observed periods with model predictions could not help distinguishing whether the star is in the PMS or post-MS evolutionary phase of its evolution, even if the solution is consistent with the distance determined by Hipparcos. The latter, if correct, implies that V351 Ori is much closer than the Orion complex.

Owing to the short temporal coverage of our previous observations, the published periodicities are subject to the problem of one cycle per day aliasing, rendering their ḏetermination rather uncertain. With the aim of overcoming this observational limitation, we have organized a multisite campaign on V351 Ori involving telescopes located in five different countries. In this paper, we present the results of this campaign. The organization of the paper is the following: in Section 2 we present details on observations and data reduction; the frequency analysis is reported in Section 3 and results are compared with those obtained in other works and with the predictions of pulsating models in Section 4; the conclusions close the paper.

## 2 Observations and data reduction

Six observatories were involved in the observational campaign. They are listed in Table 1 together with the telescopes and instruments used. The observations have been supplemented by data obtained at the SAAO during the same period (see BKW). Two comparison stars (HD 38248 and HD 290826) have been used during the observations, depending on the constraints set by the different instruments (e.g by the field of view of the CCD etc.). Some information on these objects is given in Table 2, along with the same data for V351 Ori. For completeness, the last column of Table 1 indicates the comparison star used at each telescope.

The multisite observations on V351 Ori were performed in the Johnson $`V`$ filter over around 2 years. We can divide the campaigns in: 1) winter 2001 (Sept. 2001-Feb. 2002); 2) winter 2002 (Oct. - Dec. 2002). In total, we gathered about 180 hour of observations over 29 nights. The detailed log of the observations is given in Table 3.

### 2.1 Details on single site observations

The observations obtained at SAAO are described in BKW. Considering the BAO and Loiano observations, in both cases the instrument was a three channel photometer (see www.bao.ac.cn for the BAO instrument and www.na.astro.it/$``$silvotti/TTCP.html for the TTCP in Loiano). These instruments allow the simultaneous observation of the variable and comparison stars, as well as the sky. Thus, the time sampling can be very short: for V351 Ori we adopted an exposure time of 10 s.

Observations in SPM were carried out using the six channel spectrophotometer (see 132.248.3.38/ Instruments/danes/photomdan.html) which can simultaneously operate the $`uvby`$ filters. For V351 Ori, we used only the $`y`$-band data after the proper conversion into the Johnson $`V`$ band.

Concerning CCD observations, below we report the data acquisition and reduction in detail.

### 2.2 Roque de los Muchachos Observations

The observations with the 1.0m JKT (La Palma, Canary Islands, Spain) were carried out during November-December 2001 using the JAG CCD camera equipped with the 2048x2048 pix2 SITe2 detector, providing a 11$`\times `$11 arcmin2 field of view and a pixel scale of 0.33 arcsec/pixel. A Johnson $`V`$-band filter was used during all observations with exposure times between 10 and 20 seconds (depending on airmass and seeing). When combined with the read-out time the sampling rate is decreased to about one photometric point every two and half minutes. Standard data reduction procedures were performed with IRAF. Circular aperture photometry was carried out with DAOPHOT under IRAF to find instrumental magnitudes of the simultaneously observed target and comparison stars. The adopted aperture radius was 70 pixels (23.1 arcsec) for the night of November 29, and 25 pixels (8.25 arcsec) for the nights of November 29 and December 3. Different aperture radii were used to account for the different seeing conditions.

### 2.3 Kitt Peak Observations

$`V`$-band observations were obtained with the SARA 0.9m automated telescope at Kitt Peak in Arizona. A little more than 42 hours of differential time-series photometry was obtained on nine nights between October 2001 and November 2002 by both on-site and remote-access observers. Images were collected using an Apogee AP7p camera with a back-illuminated SITe SIA 502AB 512x512 pixel CCD. The pixels are 24 microns square, corresponding to 0.73” at the telescope focal plane scale. Read out noise and gain for the camera are about 12.2 electrons (rms) and 6.1 electrons/ADU, respectively.

Integration times were typically 10-30 seconds (depending upon sky conditions and airmass), plus $``$ 5 seconds readout and idle time between frames. Sky flats, dark and bias exposures were taken every night. All data was calibrated and reduced using standard IRAF routines.

### 2.4 Teide Observations

Observations at Teide Observatory were carried out with the 1.0m OGS Telescope equipped with a 1024x1024 pixel CCD. The pixel scale was 0.32$`^{^{\prime \prime }}`$ for a total field of view of 5.5$`^{^{}}`$x5.5$`^{^{}}`$. Read noise and gain for the camera are about 5.4 electrons (rms) and 2 electrons/ADU, respectively. Typical integration times were 10-15 sec, depending upon sky conditions. Sky flats, dark and bias exposures were taken every night. All data was calibrated and reduced using standard IRAF routines.

In order to homogenize the differential photometry obtained at telescopes which used different comparison stars and to overcome the consequent difficulty in obtaining a secure zero point stability, we have decided to subtract from the data gathered during each night at each site their mean values. This procedure will limit our analysis to frequencies longer than 1 c/d. However, since the longest time series obtained at a single site lasts around 7 h, all the frequencies lower than $``$ 3.3 c/d have to be considered suspect. However, this is not a serious drawback, since the range of frequencies of interest for $`\delta `$ Scuti pulsations generally falls beyond this value.

In some case observations made at different observatories overlapped in time, in particular the run at SAAO and ING and/or Loiano. In this event, the data were averaged together. In order to compare the quality of the data from different sites taken at the same time, we show in Fig. 1 the observations collected at BAO, SAAO and ING during the night of November 29 (HJD=2452243).

In order to have an approximately similar sampling and a similar weight for each dataset in the Fourier analysis, we had to bin the photometry obtained at BAO, Loiano, OGS and SARA. As a consequence, the total number of data points has been reduced by a factor $``$4. In this way, we have the significant advantage of using much less computer time for frequency extraction.

We note that, in order to improve the analysis, and in particular the phase determination, we have added to the set of observations described above the data obtained about 1 year before and described in M01. We recall that these observations were obtained at Loiano with the TTCP. These data were binned in the same manner as described above.

Finally, the light variations of V351 Ori obtained during the winter 2001 and winter 2002 campaigns are shown in Figs. 2,3.

## 3 Frequency analysis

The frequency analysis was performed using the Period98 package (Sperl 1998, www.astro.univie.ac.at/$``$dsn/), based on the Fourier transform method. For a better interpretation of the results, we have first calculated the spectral windows (SW) for the data set collected in winter 2001, winter 2002, and for the whole data set (which includes also the 2000 M01 data). The result is shown in Fig. 4. As can be inferred from the figure, the 1 c/d alias is not significantly reduced when analysing the whole dataset with respect to the individual sets of data (from $``$80% to $``$70% in amplitude). However, the main improvement obtained from analysing the whole dataset is the strong decrease of aliases between 0 and 1 c/d.

The spectral window is useful for estimating the error on the extracted frequencies. The latter is approximately equal to the FWHM of the main lobe in the spectral window (see, e.g., Alvarez et al. 1998). In our case, we find $`\mathrm{\Delta }f`$ $``$0.03 c/d for the whole dataset, whereas the errors for the winter 2001 and winter 2002 campaigns are $`\pm 0.06`$ c/d and $`\pm 0.05`$ c/d, respectively.

Before extracting the pulsation frequencies from the time series, it is important to have a tool to identify the last significant periodicity in the data set. To this aim, we have adopted both the Scargle test (Scargle, 1982, see also Horne & Baliunas 1986) and the empirical criterion suggested by Breger et al. (1993) (see also Kuschnig et al. (1997)). These authors state that the signal-to-noise ratio (in amplitude) should be at least 4 in order to ensure that the extracted frequency is significant. In the following, we adopt this convention that will be referred to as “the Breger criterion”.

In order to apply the Scargle test, we have to estimate the noise spectrum over the whole set of investigated frequencies and the number of independent frequencies in the spectrum itself. To address the first point we have used the Period98 program to average the noise in boxes of 10 c/d width moving along the spectrum. The size of these boxes was chosen on the basis of the spectral window shape (see Fig. 4) so that at least part of the white noise was included. It is worth noting that, as recommended by Horne & Baliunas (1986), the noise spectrum used to investigate the significance of a certain frequency $`f_x`$ is calculated before subtracting that frequency.

The second point is more difficult to address, since there is no precise way of determining the number of independent frequencies in case of unevenly spaced data (which is our case). However, as discussed in Alvarez et al. (1998), we can use the FWHM of the main lobe in the spectral window ($`\mathrm{\Delta }f`$ $``$ 0.03 c/d for our dataset) to estimate this value. Using Eq. (20) of Scargle (1982), we find that 99% and 90% significance levels are reached when the signal (in amplitude) is 3.9 and 3.5 times the noise. Note that the Scargle Test is, in practice, equivalent to the Breger criterion111It is important to note that in Breger et al. (1993) the noise level is calculated after that all the relevant frequencies have been subtracted. To be conservative, we have adopted a noise level equal to the one used in Scargle’s test.. In order to rely only on secure frequency detection, we choose to consider only frequencies whose level of confidence is better than 99% (i.e. S/N=3.9).

Finally, we performed the Fourier analysis of the whole data set. This procedure identifies five frequencies that satisfy the above requirement. The periodograms are shown in Fig. 5, where the solid line displays the S/N=4 level, and the dotted and dashed lines represent the 99% and 90% significance levels of the Scargle test, respectively. The derived frequencies, amplitudes, phases and relative S/N are listed in Table 4; the corresponding phase diagrams are shown in Fig. 6. A comparison between the frequencies extracted from the whole data set and those obtained during the single winter 2001 and 2002 campaigns is given in Table 5. Inspection of the table reveals that frequencies $`f_1`$ to $`f_4`$ have been detected within the errors both in the winter 2001 and 2002 campaigns, whereas $`f_5`$ is present only in the winter 2002 campaign.

We note that in Fig. 5 there is another apparently significant frequency at 3.347 c/d with S/N$``$3.9. However, we do not consider it relevant to our discussion because the period of $``$ 7.2 hr is just the duration of the longest time series obtained at a single site (see Sect. 2), and the resulting frequency could be a spurious artifact. Whether or not this is the case goes beyond the interest of this paper, since its inclusion in the fit has negligible influence on the results for the other frequencies, and in any case the frequency is much shorter relative to the typical values of $`\delta `$ Scuti stars.

In order to test the Period98 results, we have independently carried out a frequency analysis using the Starlink Period package (star-www.rl.ac.uk). In particular, we used the Lomb Normalized Periodogram (LNP, Lomb 1976; Scargle 1982), noting that other algorithms, such as CLEAN (also included in Period) give results similar to those obtained with LNP.

The result of this analysis is shown in the third column of Table 6, labelled LNP (1). There is an obvious discrepancy concerning $`f_5`$, which is not detected by LNP. To check if $`f_5`$ is reliable or not, we have adopted the following procedure: 1) we analysed the time series by means of LNP; 2) once extracted the first frequency, we used Period98 (instead of the Period package) for the fit and prewithening; 3) the residuals from the prewithening have been analysed again using LNP; 4) we iterated the previous steps. The result of such a procedure is shown in the fourth column of (labelled as LNP (2)) of Table 6. This time the agreement with results from Period98 is almost perfect. As a conclusion, our test shows how using different packages could lead to somewhat different results, suggesting the need to use more than one method to analyse the data.

As a result of the above considerations we can conclude that the frequencies $`f_1`$, $`f_2`$, $`f_3`$, and $`f_4`$ are established beyond any reasonable doubt, whereas $`f_5`$ is slightly less reliable.

The frequency $`f_5`$ deserves more discussion. In fact, it is very close to $`f_3`$, the separation being only 0.06 c/d (we recall that the resolution for the whole data set in 0.03 c/d). As noted by Breger & Bischof (2002) close frequency pairs have important asteroseismological implications, and it is prudent to investigate whether they are real or not. In fact, there are many observational problems listed by Breger & Bischof (2002) that can lead to a false closely spaced double modes. In particular, amplitude variability of a single frequency can resemble the presence of a close pair. However, from the observational point of view it is very difficult to discern which hypothesis (close pair or single mode with amplitude variations) is correct. To this aim, Breger & Bischof (2002) illustrate a method which relies on the examination of amplitude and phase variation with time. However, this analysis requires a huge amount of observations while the data on V351 Ori presented here is insufficient for such a test and we simply warn the reader that $`f_5`$ might not be a real pulsation frequency.

## 4 Discussion

### 4.1 Comparison with previous work

Table 7 shows the frequencies derived in the present work, along with those found by M01 and BKW. With respect to M01, we are able to confirm within the uncertainties only 3 frequencies: $`f`$=15.49 c/d, $`f`$=11.89 c/d (this is likely the -1 c/d alias of our $`f_3`$), and $`f`$=16.27 c/d. Concerning the other three frequencies reported in M01, they are likely due the noise in the observations.

On the basis of the SAAO photometry and radial velocity measurements of V351 Ori, BKW have identified two and three pulsational frequencies respectively (see the last two column of Table 7). In particular they found 15.675 c/d and 14.335 c/d from the photometry and 15.682 c/d, 14.153 c/d and 11.877 c/d from radial velocities data (two of which are clearly the same as found in the photometry). BKW also performed a frequency analysis on the basis of line profile variations (see their Sect. 6 for details) that yielded a frequency at 1.90 c/d, considered highly uncertain and found indications for an additional frequency at about 20 c/d. Comparison with the BKW results (see Table 7) shows that we confirm the presence of the two frequencies at 15.675/15.682 c/d and 14.335/14.153 c/d (our $`f_1`$ and $`f_2`$). The third feature at 11.877 c/d could be interpreted as the $``$1 c/d alias of our $`f_3`$ within 2 $`\sigma `$ or of $`f_5`$ within 1 $`\sigma `$. Alternatively, since our 1 c/d alias is not negligible, our frequency $`f3`$ could be the alias and the frequency at 11.877 c/d the correct one. An additional hypothesis is that this frequency does not correspond neither to $`f_3`$ nor to $`f_5`$, i.e. it is an additional pulsational frequency which shows up only with spectroscopic measurements. Concerning the remaining two periodicities reported by BKW, we do not find any power at 20 c/d and we cannot say much on the frequency at 1.90 c/d, because of our reduction technique.

### 4.2 Theoretical constraints

Following the same approach adopted in M01, we have tried to reproduce the five $`\delta `$ Scuti-like frequencies derived in the previous section using linear non-adiabatic models of PMS $`\delta `$ Scuti stars. For each selected mass, the model luminosity and effective temperature were chosen from PMS evolutionary calculations and the constraint is provided by the observed pulsation period(s). Note that the $`f_1`$ and $`f_4`$ frequencies cannot be reproduced simultaneously since their separation is too small and comparable to the theoretical uncertainty due to the model mass zoning. The same holds for $`f_3`$ and $`f_5`$, but, as pointed out in the previous section, the latter frequency is slightly less reliable.

On this basis, we have tried to reproduce only the three frequencies with the highest amplitude, $`f_1`$, $`f_2`$, and $`f_3`$. The results show that two radial pulsation models can be found that simultaneously pulsate with frequencies $`f_1`$ and $`f_3`$, corresponding to a mixed oscillation either in the first and second overtone or in the second and third overtone. We caution that at present we cannot indicate which of the models is the real one. The stellar parameters of the two best fit pulsation models are given in Table 8 and the position in the HR diagram is displayed in Fig. 7 (filled circles). The dashed box indicates the range in luminosity and effective temperature corresponding to the empirical estimates available in the literature. The minimum luminosity comes from the Hipparcos lower limit to the distance of 210 pc (van den Ancker et al. 1998), whereas the upper limit assumes that V351 Ori is located in the Orion star forming region. The width of the box corresponds to an uncertainty of $`\pm `$0.01 dex around the mean value of $`\mathrm{log}T_{eff}3.88`$ (van den Ancker et al. 1998). The position of the best fit models indicates that the distance of V351 Ori is intermediate between the lower limit set by Hipparcos and that of the Orion Nebula.

Despite the success in reproducing $`f_1`$ and $`f_3`$, our radial models fail to also account for $`f_2`$. This fact, and the inability to distinguish between $`f_1`$ and $`f_4`$, and between $`f_3`$ and $`f_5`$ indicate that V351 Ori may be pulsating in a mixture of radial and nonradial modes, in agreement with the early suggestion of BKW. The possibility that nonradial modes might be associated with the identified frequencies clearly requires more theoretical investigation. Preliminary tests based on the comparison between observed and predicted large frequency separations (by means of scaling laws) seem to suggest that the 2$`M_{}`$ best fit model can be compatible with the observations if we consider that $`f_1`$ and $`f_3`$ are $`l=0`$ modes of consecutive mode order, whereas $`f_2`$ is a $`l=1`$ mode. A similar test for the 2.3$`M_{}`$ model suggests that in this case it is much more difficult to reconcile model predictions with the observed frequencies and frequency separations. However, these results, based on scaling considerations, need to be confirmed through the computation of specific nonradial models and we plan to address this issue in detail in a forthcoming paper.

## 5 Conclusions

A total of 180 hour of observations obtained during 29 nights at seven different telescopes on the PMS $`\delta `$ Scuti star V351 Ori have been presented. The Fourier analysis of this data set confirms the multiperiodic nature of this pulsator: we have identified five frequencies of pulsation, four of which are highly significant, while the last is more uncertain even if it is significant on the basis of both the Scargle test and the Breger Criterion.

The frequencies presented in this paper differ from those reported in M01, implying that the theoretical interpretation presented there needs to be revised. The comparison with the predictions of linear nonadiabatic radial pulsation models, in the hypothesis that the mode identification is retained, suggests that the distance to V351 Ori is intermediate between the lower limit measured by Hipparcos and the distance to the Orion Nebula. However, the fact that radial pulsation models are not able to reproduce simultaneously all the observed frequencies strongly suggests that nonradial modes might also be present. We defer a complete analysis of the pulsation modes to a forthcoming study.

###### Acknowledgements.

We are indebted to L. Balona, C. Koen and F. van Wyk for having provided us with their data in advance of publication and in particular to C. Koen for a critical reading of the manuscript.We acknowledge R. Silvotti and L. Piau for helpful discussions. We wish also to thank the Loiano staff members for their kind support during the observations. AAF is thankful to DGAPA-UNAM project IN113599 for financial support. FJGP acknowledges support from project ESO/FNU/ 43658/2001. DFMF acknowledges financial support by FCT from the “Subprograma Ciência e Tecnologia do 3o Quadro Comunitário de Apoio” and also from project POCTI/1999/FIS/34549, approved by FCT and POCTI, with funds from the European Union programme FEDER.This work was partly supported by FCT-Portugal, through project POCTI/43658/FNU/2001. TDO would like to acknowledge partial support for this project from NASA (grant NAGW5-9408) and NSF (AST-0097616).


format: text

text:
# The Transition between Nonorthogonal Polarization Modes in PSR B2016+28 at 1404 MHz

## 1 INTRODUCTION

All modern polar cap theories of pulsar radio emission are based upon the hypothesis of Radhakrishnan & Cooke (1969, hereafter RC) that the emission arises from charged particles streaming along open magnetic field lines above the pulsar’s polar cap. The observed polarization due to this particle motion is the instantaneous projection of the magnetic field lines in the emission region on the plane of the sky. If the magnetic field structure is dipolar, the position angle (PA) of the resulting linear polarization vector varies in an S-shaped pattern with pulse longitude as the pulsar’s rotation alters the field-line projection with respect to an observer’s sightline. Since the polarization vector’s orientation is determined solely by the dipolar field structure and the viewing geometry, the observed PAs are essentially frequency-independent. The astonishing ability of the RC model to explain both the PA patterns and the large fractional linear polarization that are actually observed is a testament to the model’s success and popularity for over three decades.

The 1404-MHz polarization observations of individual pulses from PSR B2016+28 by Stinebring et al. (1984, hereafter SCRWB) offer two very different interpretations of the pulsars’ viewing geometry, thereby challenging the validity of the RC model (SCRWB). The longitude-dependent PA histograms presented by SCRWB reveal two, nearly flat, PA traces separated by about $`90^{}`$ across most of the pulse with a connecting bridge between the traces near the pulse center. The two PA traces most certainly arise from the well-known orthogonal modes of polarization (e.g. Manchester, Taylor, & Huguenin 1975, hereafter MTH; Backer, Rankin, & Campbell 1976, hereafter BRC; Cordes, Rankin, & Backer 1978; Backer & Rankin 1980, hereafter BR; SCRWB), and their flat trajectories suggest a viewing geometry where the observer’s sightline makes a large angle with respect to the star’s magnetic pole. But when the SCRWB data are averaged over many rotations of the star, the resulting PAs vary in the classic S-shaped pattern with a steep slope near the pulse center (Fig. 1), suggesting that the sightline nearly intersects the magnetic pole. The steep PA sweep observed in the averaged data is embedded in the modal connecting bridge shown in the single pulse data. The rapid PA change and modal connecting bridge are not apparent in the 430-MHz polarization observations reported by BRC, as one might expect from the frequency independence of the RC model. The PAs of PSR B2016+28 observed at 1404 MHz then pose a multi-faceted dilemma. If the RC model is to be preserved, one must decide which viewing geometry is correct while also developing an explanation for the PA trace that is not attributed to geometry. Furthermore, one must determine the origin of the modal connecting bridge. Is it an artifact of polarization mode interaction or could it be an indicator of another emission mechanism operative in the pulsar magnetosphere?

Despite the conflicting interpretations generated by the 1404-MHz PAs of PSR B2016+28, its viewing geometry has generally been accepted on the basis of the frequency dependence of its pulse shape. Rankin (1983, 1986, 1993) developed a morphological classification scheme that is based primarily on the spectral evolution of pulse shapes. She categorizes PSR B2016+28 within a group of pulsars designated as “conal single” or “type-Sd” pulsars. The pulse shape of a conal single pulsar evolves from a single component at high frequency to two components at low frequency. This frequency dependence of the pulse shape is thought to result from a combination of the spreading or radius-to-frequency mapping of the conal emission and a tangential traverse of an observer’s sightline across the pulsar’s emission cone. The viewing geometry for conal single pulsars thus requires a large angle between the magnetic pole and the sightline, just as the modal PA trajectories in PSR B2016+28 suggest. The pulsar displays other key characteristics of conal single pulsars (Rankin 1983; Rankin & Ramachandran 2003), such as drifting subpulses and low linear polarization across its entire pulse owing to orthogonally polarized modes (OPM). Hankins & Rickett (1986) investigated the frequency dependence of the average profiles for many pulsars. Specifically recognizing the conflicting interpretations of the 1404-MHz PAs in PSR B2016+28, they noted that the modal PA trajectories in the pulsar are more consistent with the frequency dependence of its profile. Hankins et al. (1992) interpret the 1404-MHz PAs, and thus the pulsar’s viewing geometry, in a similar way. However, no one has adequately resolved the dilemma over the conflicting interpretations of the 1404-MHz PAs or reconciled, in detail, the 1404-MHz PAs with the viewing geometry implied by its morphological classification. Additionally, no one has offered a detailed explanation for the pulsar’s modal connecting bridge.

SCRWB suggested that the modal connecting bridge in PSR B2016+28 and, more generally, the nonorthogonality of the polarization modes may be related to the “migration” of OPM. This mode migration may be the frequency-dependent relocation, spreading, or splitting of modal radiation with pulse longitude. SCRWB also emphasized the need to clarify the nature and importance of nonorthogonal polarization modes (NPM) in the emission. While the occurrence of NPM has been frequently observed and widely documented in the literature (e.g. BR; SCRWB), theoretical investigations (e.g. Cheng & Ruderman 1979; Allen & Melrose 1982; Barnard & Arons 1986; Michel 1987; Gangadhara 1997) and statistical modeling (e.g. McKinnon & Stinebring 1998, 2000; McKinnon 2002) of the modes have concentrated on their orthogonality, but have largely ignored their observed nonorthogonality.

The first objective of this paper is to resolve the dilemma over the conflicting interpretations of the pulsar’s 1404-MHz PAs. I follow the suggestion of SCRWB and argue that the rapid change in average PA along with the modal connecting bridge in PSR B2016+28 is a longitudinally-resolved transition between NPM. In other words, the bridge is simply an artifact of NPM superposition, and not an indicator of another emission mechanism. The second objective of the paper is to account for nonorthogonal polarization modes in a statistical model of the radio emission’s polarization. In §2, I review the observed polarization properties of PSR B2016+28 by replicating the results of SCRWB, and argue that the modal PA trajectories reflect the true viewing geometry of the pulsar. I incorporate the mode nonorthogonality in a statistical model of superposed OPM (McKinnon & Stinebring 1998, hereafter MS1) in §3. I then describe the predicted effects of superposed NPM on polarization data and compare the predictions with observations. In §4, I discuss the implications of this investigation for radiative processes in pulsar magnetospheres.

## 2 OBSERVED POSITION ANGLES IN PSR B2016+28

As a point of reference for the discussion in § 3, the observational results of SCRWB are replicated in Figure 1, which shows the average pulse profile of PSR B2016+28 at 1404 MHz. The single pulse, polarization data used to construct the profile were recorded by SCRWB with the 305-m Arecibo radio telescope. The solid line in the bottom panel of the figure shows how the average PA varies across the pulse. The PAs of the pulsar’s nearly-orthogonal polarization modes at each pulse longitude are shown by the triangles in the PA plot, and were estimated as follows.

A PA histogram, similar to those shown in Figure 5, was computed at each pulse longitude of PSR B2016+28 using the single pulse, polarization observations of SCRWB. Each histogram contained 50 equally-spaced bins over the $`180^{}`$ range in possible PA values, giving a histogram bin resolution of $`3.6^{}`$. Only data samples having a linear polarization that exceeded the off-pulse noise by a factor of five were used to construct the histograms. The PA of the most frequently occurring, or dominant, polarization mode was found by locating the peak in the histogram, and using the peak bin along with 12 bins on either side of it to calculate a weighted mean angle, $`\psi `$. From the statistics of directional data described in Mardia (1972), the weighted mean angle was found using

$$C=\frac{1}{N}\underset{i=1}{\overset{25}{}}f_i\mathrm{cos}(2\psi _i),$$
(1)
$$S=\frac{1}{N}\underset{i=1}{\overset{25}{}}f_i\mathrm{sin}(2\psi _i),$$
(2)
$$\psi =\frac{1}{2}\mathrm{arctan}\left(\frac{S}{C}\right),$$
(3)
where $`\psi _i`$ is the midpoint of the PA bin, $`f_i`$ is the number of data samples in the bin, and $`N=f_i`$ is the total number of data samples in the 25 PA bins. The open triangles in Figure 1 show the weighted mean angles for the dominant mode at each pulse longitude. The remaining 25 bins in the PA histogram were used to calculate the mean angle of the weak polarization mode (filled triangles in the figure) with the same equations. The resulting PA trajectories of the polarization modes are very similar to those reported in the original observations (see Fig. 31 of SCRWB).

Of the 2700 pulses recorded in the observation, the number of data samples exceeding the detection threshold on linear polarization was as high as 1250, depending on pulse longitude. Data points are not shown in the bottom panel of Figure 1 for any pulse longitude having less than 100 samples in its PA histogram. These locations occur at the pulse edges where the instantaneous signal-to-noise ratio (S/N) in linear polarization is low. Formal errors in $`\psi `$ are estimated to range from $`6^{}`$ to $`13^{}`$, and have been omitted from the figure for clarity.

The primary advantage of calculating modal PAs from equations 1, 2, and 3 is to bring to bear the full weight of the data in making quantitative estimates of modal PAs and their differences, thus helping to identify any deviations from mode orthogonality. A minor disadvantage of this particular method is to introduce a slight bias that skews the difference in modal PAs towards $`90^{}`$ when the PA is nearly uniformly distributed and the linear polarization is very low (e.g. on the pulse edges and near a pulse longitude of $`1^{}`$ in Fig. 1).

The PAs in the bottom panel of Figure 1 clearly illustrate the two interpretations that could be made of the pulsar’s viewing geometry. The average angle changes rapidly with pulse longitude near the pulse center, and its total excursion across the pulse approaches $`180^{}`$. Both properties are indicative of a viewing geometry where the observer’s sightline passes near the star’s magnetic pole. The trajectories followed by the polarization modes, however, are relatively flat with a total excursion in PA of something much less than $`180^{}`$, indicating that the observer’s sightline is far from the magnetic pole.

Many features of Figure 1 suggest that the flat PA trajectories of the polarization modes reflect the pulsar’s true viewing geometry. The average PA generally follows the dominant polarization mode (open triangles in Fig. 1) as one would expect from OPM superimposed upon a flat PA trajectory (BRC; SCRWB; MS1). One polarization mode dominates the leading edge of the pulse, while the other dominates the trailing edge. Assuming that the modes occur simultaneously (MS1; McKinnon & Stinebring 2000, hereafter MS2), the low average linear polarization and the frequent occurrence of both modes across the entire pulse imply that the modes have comparable polarized intensities on average (see, also, Rankin 1983; Rankin & Ramachandran 2003). Additionally, the rapid change in average angle near the pulse center isn’t quite as smooth as one would typically expect from the RC model. The rapid PA change joins the modal PA trajectories, and its magnitude is approximately $`90^{}`$, suggesting that the PA change is related to an OPM transition.

Additional features of Figure 1 suggest that NPM may be complicating our interpretation of the PA pattern in PSR B2016+28. The differences between the modal PAs are almost, but not precisely, $`90^{}`$ across most of the pulse. The mode nonorthogonality does not appear to vary greatly or rapidly across the pulse. While the change in the average angle at the pulse center is rapid, it is not discontinuous as would be expected for a purely orthogonal transition, and so may indicate the presence of NPM (Cheng & Ruderman 1979; SCRWB; Michel 1987).

## 3 STATISTICAL MODEL OF NONORTHOGONAL MODES

A statistical model for the polarization of pulsar radio emission was presented in MS1 and further developed in MS2 and McKinnon (2002). The model is based upon the observation that the highly polarized, orthogonal modes occur simultaneously. To incorporate the nonorthogonality of the modes in the model, the analysis which follows concentrates on the emission’s linear polarization and neglects its circular polarization because NPM are most evident in the observed PA of the linear polarization vector and the circular polarization generally tends to be a small part of the emission’s polarization.

As in MS1, let $`X_1`$ and $`X_2`$ be random variables representing the linearly polarized intensities of the primary and secondary polarization modes, respectively. Without loss of generality, the PA of the primary mode can be set to $`0^{}`$, which means that the polarization vector of the primary mode points in the direction of increasing Stokes parameter Q in the Q-U plane of the Poincaré sphere (Fig. 2 ). The polarization vector of the secondary mode generally points in the direction opposite that of the primary mode vector, but makes an angle $`\theta `$ with respect to the Q-axis to account for the nonorthogonality (see Fig. 2). If $`X_\mathrm{N}`$ is a random variable representing the instrumental noise, the observed Stokes parameters Q and U of the combined emission at an arbitrary pulse longitude are

$$\mathrm{Q}=X_1X_2\mathrm{cos}\theta +X_{\mathrm{N},\mathrm{Q}},$$
(4)
$$\mathrm{U}=X_2\mathrm{sin}\theta +X_{\mathrm{N},\mathrm{U}}.$$
(5)
Assuming that $`\theta `$ is constant and thus independent of $`X_2`$, the PA calculated from the average values of Q and U is

$$\psi =\frac{1}{2}\mathrm{arctan}\left(\frac{\mathrm{U}}{\mathrm{Q}}\right)=\frac{1}{2}\mathrm{arctan}\left(\frac{\mathrm{sin}\theta }{M\mathrm{cos}\theta }\right),$$
(6)
where $`M=\mu _1/\mu _2`$ is the ratio of the mean values of the modes’ polarized intensities, $`\mu _1=X_1`$ and $`\mu _2=X_2`$. Equation 6 shows that the average PA resulting from the superposition of NPM produces an angle intermediate to the PAs of the two modes. For example, when $`M=1`$ and the nonorthogonal modes have PAs of $`0^{}`$ and $`85^{}`$, the deviation from orthogonality is $`\theta =10^{}`$ and the average angle is $`\psi =42.5^{}`$, which is equidistant from the orientations of the two modes.

If the modes are truly orthogonal ($`\theta =0`$), the vectors representing their linear polarization are antiparallel in the Q-U plane. The polarization vector resulting from the superposition of OPM is aligned with the PA of the mode having the larger polarized intensity (Fig. 2, example 1). The MS1 statistical model attributes the observed switching between modes to the random fluctuations in their polarized intensities. Within the context of the model, the only way to produce PAs intermediate to those of the orthogonal modes is through the influence of the instrumental noise. PA histograms produced by this model generally contain two peaks, symmetric about their centroids and separated by $`90^{}`$, with a uniform plateau arising from the instrumental noise (see Fig. 2 of MS1).

When the modes are not orthogonal, their polarization vectors are not antiparallel in the Q-U plane. As can be seen in examples 2 through 4 of Figure 2, the orientation of the vector resulting from NPM superposition is intermediate to the orientations of the mode vectors, and depends upon the modes’ polarized intensities and the degree of nonorthogonality ($`\theta `$ in Fig. 2). Also, the superposition of NPM depolarizes the observed emission (i.e. the amplitude of the resultant polarization vector is less than that of the polarization vector for the dominant mode), but not nearly as efficiently as the superposition of OPM.

These signatures of the superposed NPM model should be evident in PA histograms and in the joint probability density, or scatter plots, of PA and linear polarization (L). The PA histograms should be similar to Figure 3, which was constructed from a numerical simulation of equations 4 and 5. The separation of the peaks in the histogram is less than $`90^{}`$ because the modes are not orthogonal. Since PAs intermediate to the orientations of the mode vectors can result from the superposition of NPM, an overabundance of data samples occurs between the peaks, and a paucity of data points occurs outside them. The intermediate values of the resultant PAs also lead to slight asymmetries in the peaks about their centroids. For the L-PA scatter plots, the NPM model predicts that data samples with significant L should occur between the mode peaks.

Figure 4 shows L-PA scatter plots for four pulsars observed by SCRWB. In each plot, the dominant polarization mode is located where the concentration of data points is greatest. The observed Stokes parameters Q and U of each data point have been rotated so that the dominant mode occurs at a PA of $`\psi =45^{}`$ for each pulsar. As predicted by the superposed NPM model, many data samples with highly significant linear polarization ($`\mathrm{L}>5\sigma _N`$) are located between the mode peaks ($`45^{}<\psi <45^{}`$). This is particularly true of PSR B0950+08 and PSR B1929+10 where the modes are clearly defined. The modes are not as well defined in the scatter plots for PSR B2016+28 and PSR B2020+28 because the modes occur with similar frequency. As discussed in MS1, this implies the mode polarization amplitudes are similar on average (i.e. $`M1`$), so that their simultaneous interaction leads to significant depolarization of the intrinsic radiation and causes the observed polarization to be dominated by the randomly polarized, on-pulse, instrumental noise (the broad pedestal in the scatter plots).

Figure 5 shows the PA histograms constructed from the scatter plots in Figure 4. As mentioned in § 2, only data points with a S/N in L greater than five (the dashed horizontal line in each plot of Fig. 4) were used to construct the histograms. The separation between the two peaks in each histogram is less than $`90^{}`$, and of the four examples shown in the figure, the deviation from mode orthogonality is greatest for PSR B2016+28. The four histograms show other similarities in addition to the nonorthogonality of the modes. In each histogram, more data points fall between the peaks than outside them, and the individual peaks are slightly asymmetric about their centroids. The PA histograms bear a remarkable resemblance to the histogram generated with the numerical simulation of the superposed NPM model (Fig. 3).

Turning now to the PA calculated from the averaged Stokes parameters, consider a hypothetical situation where one had apriori knowledge of a pulsar’s viewing geometry. This being the case and barring the limitations of S/N, one could easily remove the systematic effect of the RC PA sweep from the data. Having done this, the variation in PA across the pulse attributable to the polarization modes would be described by equation 6, where $`M`$ and $`\theta `$ are now dependent on pulse longitude. The primary mode would occur at pulse longitudes having $`M>\mathrm{cos}\theta `$, the secondary mode would occur at longitudes with $`M<\mathrm{cos}\theta `$, and a mode transition would occur at longitudes where $`M=\mathrm{cos}\theta `$. If the deviation from orthogonality is small and constant ($`\theta =\theta _o<<1`$) over the range of pulse longitude where the mode transition occurs and if the modal ratio, $`M`$, varies linearly with longitude (e.g. $`M(\varphi )=1+a(\varphi \varphi _o)`$) over the same range, the rate the PA changes with pulse longitude near the mode transition is

$$\frac{d\psi }{d\varphi }\frac{a}{2\theta _o}.$$
(7)
Equation 7 shows that the mode transition is abrupt when the modes are orthogonal ($`\theta _o=0`$), but gradual when the modes are not orthogonal ($`\theta _o0`$). The equation also suggests that a mode transition may be longitudinally resolved provided that $`\theta `$ and $`M`$ vary slowly with pulse longitude, which is likely the case for PSR B2016+28 for reasons cited in § 2.

The PA change with pulse longitude that is described by equation 6 is shown in Figure 6 for different values of $`\theta _o`$, with $`M`$ varying linearly with longitude. The change in average PA predicted by the statistical model (Fig. 6) is qualitatively consistent with what is observed in PSR B2016+28 (cf. bottom panel of Fig. 1). Equation 6 and Figure 6 also illustrate the point that a PA transition between purely orthogonal modes will be discontinuous and instantaneous, whereas nonorthogonal modes will produce a smooth or gradual PA transition (Cheng & Ruderman 1979; SCRWB; Michel 1987).

If, as in SCRWB, one were to convert the PA histograms at each pulse longitude (e.g. Fig. 3) to grey-scale histograms and overlay them on the averaged PA data (e.g. Fig. 6), the average angle would be embedded in the pedestal that appears between the peaks in each histogram because the PA resulting from the superposition of NPM is intermediate to the mode PAs (Eq. 6). In overlaying all the histograms on the average data, the individual pedestals combine to form a much broader pedestal that spans a large range in pulse longitude. This broad pedestal is the modal connecting bridge in PSR B2016+28 discovered by SCRWB (see Fig. 31 of SCRWB).

The preceding point is also illustrated in Figure 4 where actual average angles are shown by the dashed vertical lines in each L-PA scatter plot. In each case, the average angle is not precisely aligned with the dominant mode PA ($`\psi =45^{}`$). The difference between the dominant mode PA and the average PA, $`\mathrm{\Delta }\psi `$, is greatest for PSR B2016+28 and PSR B2020+28 where the modes occur with similar frequency and the mode nonorthogonality is significant. This is to be expected from the model because equation 6 predicts that $`\mathrm{\Delta }\psi `$ is large for $`M1`$ and $`\theta >0`$.

## 4 DISCUSSION

The statistical model’s ability to simulate the rich variety of polarization properties observed in pulsar radio emission provides strong support for the model’s underlying assumption of superposed polarization modes. Observations of PSR B2020+28 support the model’s prediction that the emission will be depolarized on all timescales and that the modes will occur with nearly equal frequency when the polarized intensities of the orthogonal modes are comparable (MS1; MS2). Additionally, the model adequately reproduces observed histograms of total intensity, linear polarization, fractional linear polarization, and PA in PSR B2020+28 (MS1; MS2). Since the polarization modes are superposed and elliptically polarized in general, the model also predicts that histograms of fractional circular polarization should be unimodal, in accord with the observations (McKinnon 2002). If the modes occurred separately, the histograms of fractional circular polarization would be bimodal, which is not observed. In this paper, the model resolves the longstanding dilemma over the interpretation of the 1404-MHz PAs in PSR B2016+28, accounts for its modal connecting bridge, and explains many subtle features of PA histograms by incorporating the nonorthogonality of the modes.

Obvious questions that arise from an investigation of NPM are: “what is their origin” and “how does the statistical model comport with the explanation of NPM generation”? The task of explaining the origin of OPM is daunting enough by itself, and little attention has been devoted to NPM’s origin because it has been regarded as a minor or subtle aspect of OPM generation. OPM may be intrinsic to the emission mechanism (e.g. Gangadhara 1997; Cheng & Ruderman 1979) or produced at widely separated regions which happen to have magnetic field structure that is orthogonal when projected on the plane of the sky (MTH; Michel 1987). However, the most widely accepted explanation for OPM appeals to propagation effects in the magnetized plasma above the pulsar polar cap (e.g. Allen & Melrose 1982; Barnard & Arons 1986). The radio emission can be spatially separated into two modes of orthogonal polarization due to the difference in their indices of refraction. If the refraction is strong, the emission in one mode that is generated at one location within the magnetosphere may overlap the emission from the other mode that is generated at a different location. In this particular scenario which is entirely consistent with SCRWB’s suggestion that NPM arise from OPM migration, the polarizations of the overlapping modes are not orthogonal because the modes were not generated at the same location. The observed polarization resulting from the overlap of the two nonorthogonal modes is precisely what is illustrated in Figure 2 and described by equations 4 and 5. Therefore, the statistical model is consistent with the proposition that OPM and NPM arise from propagation effects in the pulsar magnetosphere. More generally, the model is consistent with any mode production mechanism that involves the incoherent superposition of the modes. The exact details of how OPM migration actually occurs in PSR B2016+28 are unknown, but may be related to the pulsar’s drifting subpulses (Drake & Craft 1968; Backer 1973) and the occurrence of OPM within them, as in other conal single pulsars such as PSR B0943+10 (Deshpande & Rankin 2001) and PSR B0809+74 (Ramachandran et al. 2002).

In an admirable first attempt to explain NPM, SCRWB modeled the emission’s polarization as the sum of a nonmodal component and OPM. While SCRWB’s model cannot be ruled out, it introduces an additional level of complexity for radiative processes in pulsar magnetospheres by requiring the generation of both a nonmodal emission component and OPM. The statistical model, by contrast, portrays a much simpler situation where radiative processes are only required to produce two polarization modes, either at emission or through propagation, with slightly nonorthogonal polarizations. Furthermore, the model leads to a simple explanation for the modal connecting bridge in PSR B2016+28. Any other explanation for the bridge would likely resort to rather ad hoc assumptions about its origin, placing additional demands on theoretical explanations of radiative processes in the magnetosphere as a result.

A number of observations have posed serious challenges to the validity of the RC model, but the model has generally withstood the test of the observations. The polarization observations of the drifting subpulses in PSR B0809+74 by Taylor et al. (1971) indicated that PA variations were synchronized with the subpulses, and not fixed to the rotation of the star as the RC model predicted. Only recently have Ramachandran et al. (2002) shown that the PA variations within the subpulses are due to OPM transitions, and that the individual modes follow PA trajectories that are consistent with the RC model. Initial observations of OPM (e.g. MTH) also challenged the RC model because it could not account for the production of the modes. However, BRC noted from their single pulse, polarization observations of PSR B2016+28 that the RC model was relevant because the modal PA traces followed parallel trajectories that were consistent with the model. SCRWB pointed out that the 1404-MHz PAs of PSR B2016+28 could be interpreted in the context of the RC model with two, very different viewing geometries. In this paper, the PA interpretation dilemma is resolved and the RC model is thus preserved by accounting for the nonorthogonality of the pulsar’s polarization modes.

The occurrence of NPM has important implications for the interpretation of pulsar polarization and viewing geometry. As shown in SCRWB and this paper, the occurrence of NPM in PSR B2016+28 so disrupts the PA sweep that two, very different interpretations can be made of its viewing geometry. Everett & Weisberg (2001) used some of the most sensitive polarization observations recorded to date in an attempt to tightly constrain the viewing geometry of a number of nearby pulsars. Interestingly, the best fits of their data to the RC model were obtained by completely discarding data recorded on the pulse and relying on data of lower S/N recorded off the pulse. The viewing geometries they determined for PSR B0950+08 and PSR B1929+10 were very different from results obtained in previous work. The fact that both pulsars exhibit NPM (Fig. 5; BR; SCRWB) may explain why data recorded on their pulses deviated significantly from the RC model. Deviations from the RC model have traditionally been attributed to multi-polar magnetic field structure, implying that the emission region is close to the stellar surface. But if the deviations from the RC model are caused by NPM, as suggested here, the emission may originate at a higher location in the magnetosphere where the star’s dipole component dominates the multi-pole components of its magnetic field. Considering the substantial impact of NPM on the PA sweep in PSR B2016+28, propagation effects may also play a more significant role in determining the observed polarization than other mechanisms, such as relativistic aberration (Blaskiewicz, Cordes, & Wasserman 1991) and current flow above the polar cap (Hibschman & Arons 2001).

In summary, I have incorporated the nonorthogonality of the polarization modes in a statistical model of the polarization of pulsar radio emission. I used the model to show that the rapid change in average PA near the pulse center of PSR B2016+28 at 1404 MHz may be a longitudinally-resolved transition betweeen modes of nonorthogonal polarization, and argued that the PA trajectories of the polarization modes reflect the true viewing geometry of the pulsar. This interpretation preserves the RC model and avoids the unnecessary and unpalatable complication of invoking an additional radiation emission mechanism to explain the pulsar’s modal connecting bridge. The occurrence of NPM in other pulsars is bound to cause their PA sweeps to deviate from the RC model, and although largely ignored in the past, NPM should be taken into account when determining pulsar viewing geometries in future work. The continuing ability of the statistical model to simulate the rich variety in the observed polarization properties of pulsars lends support to the model’s applicability and its fundamental, underlying assumption that the polarization modes occur simultaneously.

I thank Dan Stinebring for providing the data used in this analysis, and Joanna Rankin for constructive comments on the manuscript. The Arecibo Observatory is operated by Cornell University under cooperative agreement with the National Science Foundation.


text:
# Making refactoring decisions in large-scale Java systems: an empirical stance

### Abstract

Decisions on which classes to refactor are fraught with difficulty. The problem of identifying candidate classes becomes acute when confronted with large systems comprising hundreds or thousands of classes. In this paper, we describe a metric by which key classes, and hence candidates for refactoring, can be identified. Measures quantifying the usage of two forms of coupling, inheritance and aggregation, together with two other class features (number of methods and attributes) were extracted from the source code of three large Java systems. Our research shows that metrics from other research domains can be adapted to the software engineering process. Substantial differences were found between each of the systems in terms of the key classes identified and hence opportunities for refactoring those classes varied between those systems.

## 1 Introduction

The term refactoring refers to a technique for improving code quality by making changes to the internal structure of the software without changing its external behaviour. Refactoring can be used to improve software design by, for example, moving code between classes, extracting code into new methods or classes or altering the position of classes in an inheritance hierarchy. Refactoring therefore leads the programmer to work more deeply on understanding what the code does and is thus an aid to maintenance and reuse .

The potential benefits of carrying out refactoring are reduced duplication of code, improved readability, faster development and fewer bugs. Whilst there has been considerable interest in refactoring principles, relatively little research seems to have focused on the identification of candidate classes for the refactoring process itself.

Extreme programming (XP) practitioners use refactoring to remove duplication whenever program features are added. Beck’s advice is that programmers should not “refactor on speculation …refactor when the system asks you to” . Whilst this advice may be appropriate for XP projects, refactoring of key classes may also be required when developing libraries for other programmers to use or when developing software in large teams.

We conjecture that there are certain classes in every OO system which are of such importance (in terms of the features they possess) that they should be a priority for refactoring effort. In this paper we show how these classes may be identified through two countable measures (number of methods and number of attributes) and through a Web-based graph metric which identifies the extent of each class’ coupling via inheritance and aggregation. The potential gain (PG) metric was first used as a means of identifying key web pages in a Web search and navigation system. The potential gain has been used here to rank classes in terms of their coupling patterns. By ranking classes according to all these features, candidate key classes can be identified.

## 2 Motivation and related work

The motivation for the research in this paper stems from a number of sources. Firstly, developers will inevitably want to quickly identify which classes are key to the system as a whole. This may be because they wish to avoid applying refactoring effort to such classes, and the associated problems of re-testing. Alternatively, it may be that those classes are the ones which need to be maintained most rigorously because they are key classes and exert considerable influence in the system as a whole. The research in this paper allows those classes to be identified more easily.

Another motivation for the work in this paper stems from a lack of research into refactoring decision criteria. The basis on which a particular class (or refactoring) should be chosen is not well understood or documented. We view it as a fruitful and interesting research topic. While there has been plenty of evidence for the benefits of refactoring, exactly what to refactor (and when) is still an open issue. Availability of tool support for refactoring (although improving significantly) seems dwarfed by current progress and thinking in the area. Our research is an attempt to highlight the criteria by which refactoring effort can be initially chanelled.

A further motivation for the research is to identify whether Java classes exhibit similar properties to their C++ counterparts. This builds on earlier work , where high-level class metrics were collected from the Unified Modelling Language (UML)  documentation of five large C++ systems, two of which were libraries. Two conjectures were investigated to determine features of key classes in a system and to investigate any differences between library-based systems and other systems in terms of coupling. Key classes in the three application-based systems tended to contain significant amounts of aggregation, large numbers of methods, attributes and associations, but little inheritance. No consistent, identifiable key features could be found in the two library-based systems; both showed a distinct lack of any form of coupling (including inheritance) other than through the C++ friend facility.

In terms of seminal refactoring literature, the PhD. work of Opdyke  describes a number of software refactorings. This thesis spawned a large amount of research in the subject. Opdyke and Johnson  describe a study in which they illustrate how to create abstract superclasses from other classes by refactoring. They decompose the operation into a set of refactoring steps, and provide examples. They also discuss a technique that can automate these steps making the process of refactoring much easier. In Johnson and Opdyke , some common refactorings based on aggregation are reported, including how to convert from inheritance to aggregation and how to reorganize an aggregate/component hierarchy. They also describe how to refine aggregations by moving variables and functions between aggregate and component classes, and how to move variables and functions within inheritance hierarchies. A text by Fowler  describes seventy-two types of refactoring and illustrates each type with examples and UML notation.

Recent empirical work in the refactoring area and its automation is found in , where 14 000 lines of code were transformed automatically where they would otherwise have had to be coded by hand. In Najjar et al. , the opportunities, benefits and problems of refactoring class constructors across a sample of classes from five Java systems were investigated. The refactoring used, replacing multiple constructors with creation methods was applied to each of a set of classes containing three or more constructors.

## 3 Empirical Evaluation

We address the empirical evaluation through a single conjecture, related to key classes as follows:

* In the type of systems investigated, a certain number of classes have higher numbers of methods, attributes, aggregation relationships and subclasses than other classes in the same type of system. These classes will, typically, be found towards the root of an inheritance hierarchy so that subclasses may take advantage of the large functionality and features they offer.

This is what we believe characterises a key class. We accept that there are many other forms of coupling which are equally worthy of investigation, and hence we do not claim our criteria to be definitive. In Section 5, we justify our choice of criteria by reference to several core refactorings.

### 3.1 Potential Gain

The potential gain metric was originally developed to aid the selection of pages allowing for future navigation of Web sites. By considering the number and length of all possible paths from a given node, an estimate of the utility or connectivity of a node could be established. If the density of the neighbourhood of some node in a graph is high, then many nodes are connected via short paths and the PG of that node is also higher. In complex systems, the lengths of these paths reflect the likely influence on connected nodes or classes.

As part of this research, a system called AutoCode was developed for indexing Java source code. AutoCode works by using a custom doclet which extends the Javadoc program and allows easy access to the code structure. We used the AutoCode system to generate graphs for each of five coupling types - Inheritance, Aggregation, Interface, Parameter Type and Return Type. Previous work described features of these graphs . This paper describes new work investigating the utility of the potential gain metric for analysing key classes using the Inheritance and Aggregation graphs.

The *Potential Gain*, PG$`(c)`$, of a class, $`c`$, is defined as the sum for all lengths (or depths) of the product of the fraction of all possible coupling paths which are of length $`d`$ and a discounting function $`f(d)`$. This gives a measure of the number of objects with which an instance of the class might potentially interact. For example, in the case of inheritance, the potential gain metric reflects both the breadth and depth of a hierarchy. A formal definition of potential gain is given in appendix A.

Related to class coupling, the PG metric gives an indication of the inter-connectivity of a class with other classes. The type of connectivity can be defined in the software which extracts (in the case of this research) the aggregation and inheritance relationships. For example, a high PG value for a class at a point in the inheritance hierarchy means that the class has a relatively large number of descendants. A very low value indicates that it is a leaf node - in other words, it has no subclasses.

Although the mathematics behind the potential gain is non-trivial, computation of the values is efficient. For example, once an appropriate graph structure had been obtained, PG values were computed for all 6 000 classes in the JDK in less than five seconds on a desktop-class PC.

### 3.2 Data collected

Data was collected from three large Java systems. Those systems were chosen since they were the subject of previous research. The three systems chosen are used extensively by developers in industrial settings. In keeping with previous work on C++ systems, we investigate two applications and one library-based system. As well as identifying key classes, previous work also found distinct differences between these two types of system, together with a lack of use of inheritance in all types of system. While this latter feature is not likely to be the case in a Java setting (because of the emphasis on inheritance in Java), our belief is that Java classes will exhibit some contradicting features across the systems. The three systems investigated were:

1. JDK \- The core Java class libraries shipped with the Java Developers Kit (JDK) which provide implementations of common functions required for many programs. This contains 1 400 000 lines of code spread over 6 000 classes.
2. Tomcat is a servlet container used in the official reference implementation for Java Servlets and JavaServer Pages. The source code for Jakarta Tomcat contains 150 000 lines spread over 370 classes.
3. Ant is a Java-based build tool and behaves in a similar way to make but uses XML-based configuration files defining various tasks to be executed. The source code for Apache Ant contains 145 000 lines of code spread over 500 classes.

#### 3.2.1 Metrics collected

The following data was collected automatically, using the same software which produced values for the PG metric.

1. The Number of Methods in a class, including public, protected and private member functions.
2. The Number of Class Attributes, including public, protected and private declarations.
3. The Depth or level of a class in an inheritance hierarchy where a zero value represents the root. The metric is based on the Depth of Inheritance Tree metric of Chidamber and Kemerer .

## 4 Data Analysis

### 4.1 Summary Data

Table 1 gives the summary data for the three systems, in terms of maximum and median values for the three metrics collected. The JDK shows the largest maximum values for all three metrics. The class with 254 methods is java.awt.Component with 80 attributes and is found at level 1 in the JDK inheritance hierarchy. A Component is any object having a graphical representation and that can interact with the user. The JDK class with 329 attributes is class xalan.templates.Constants with zero methods and again, was found at depth 1. Two classes PIORB and NSORB reside at the lowest depth of the JDK library (depth 8) with 41 and 2 methods, respectively.

The Tomcat and Ant systems are comparable in terms of number of methods and attributes. The Tomcat class with 168 methods is StandardContext, with 49 attributes and is found at depth 2. The class with 57 attributes was JspC with 45 methods, found at depth 1. This class provides the shell for the jspc compiler.

The Ant class with 89 methods is Project, found at depth 1 with 33 attributes. The class, CBZip2InputStream with 47 attributes has 31 methods (see Table 1) and is a possible candidate for refactoring. Visual inspection of both classes revealed two bad smells - Primitive Obsession and Long Methods . In our analysis we have viewed basic system types, such as String as primitive.

### 4.2 Investigating Conjecture C1

Identification of key classes has significant implications for refactoring and ultimately, maintainability. If such classes contain significant functionality and are coupled to numerous other classes, then modification of those classes needs to be made with particular care. This is particularly important if they are found at shallow levels in an inheritance hierarchy (i.e., a low depth) since all subclasses need to be considered if any change is made to that class.

From the summary data and the previous discussion, it appears that in each of the three systems there are classes with large numbers of methods and attributes. Conjecture C1 attempts to clarify the extent to which combinations of all features, including aggregation, are found in the same classes. By aggregation relationships, we mean classes which either use a large number of other classes or are used by a similarly large number of other classes (or both). Hereafter, we will call these two types of aggregation normal aggregation and reverse aggregation, respectively.

#### 4.2.1 The JDK library

Table 2 shows, for the JDK system, the numbers of methods and attributes and the depths of inheritance for the top fifteen classes when ranked in descending order according to their reverse aggregation PG value. A high reverse aggregation PG value implies that a lot of classes use the class in question. Classes such as PageAttributes.MediaType, Color and Character.UnicodeBlock are used in constant (static final) declarations and are often self-referencing.

Two of the classes identified in Table 2 have zero attributes and two have zero methods. The maximum depth in the inheritance hierarchy was 3 for two classes - Vector and javax.print.attribute.standard.MediaSizeName. The presence of Hashtable and Vector as commonly used objects suggests that the Collections framework introduced in JDK 1.2 has not been fully adopted within the JDK.

The top fifteen classes were then ranked in descending order according to their normal aggregation PG value. These fifteen classes were then compared with the reverse aggregation classes in Table 2. Eight classes were found to coincide with the previous fifteen found. Table 3 names these classes, together with the position they occupied in Table 2.

The classes in Table 3 are mostly self-referencing classes with many static references. For example, PageAttributes.MediaType contains a self-reference for each paper format (e.g. A4, US Letter, etc.). Such classes induce high PG values. A similar effect was observed with Web-search metrics. Lempel and Moran showed that metrics such as HITS  and PageRank  can be influenced by a phenomenon called the Tightly Knit Community (TKC) Effect . The TKC effect has a negative influence on search results when small groups of pages are heavily interlinked. In the analysis of program code, the TKC effect is useful for identifying classes or networks of classes which are tightly bound and where functionality is impossible to extract. This is made possible by considering more than one metric in the analysis.

A high normal aggregation PG value would imply that an instance of this class will make many references to objects which, in turn, make many references to other objects. Modifying or refactoring class with high normal aggregation PG values require careful thought and preparation, since many other classes may be affected by such a change.

By considering those classes which score highly on both metrics, we can eliminate those classes which are self-referencing. The results of eliminating such classes leave String, Object, Hashtable, Vector and Class clearly identifiable as extensively used, key classes.

Table 4 shows the top fifteen classes when ranked on descending inheritance PG values. It is interesting to note that only one class from Table 3 appears in the top fifteen JDK classes from Table 4. This class was java.lang.Object, as might be expected; a high inheritance PG value implies that a class has many subclasses.

We note that the JDK class java.awt.Component with 254 methods (Table 1) was ranked 25th in terms of reverse aggregation, 26th by aggregation and 5th by inheritance. This places the class in the top 1% for all metrics. We would therefore consider this a key class. Additionally, it may be considered a Large Class  and hence a candidate for refactoring via Extract Class or Extract Subclass because of this “bad smell”.

Finally, it is also noticeable that six classes with the highest reverse aggregation PG values have three or more constructors. This would make those classes eligible for refactoring by replacing multiple constructors with creation methods as suggested by Kerievsky and empirically investigated by Najjar et al. .

Considering conjecture C1, the classes most often used in aggregation relationships do not necessarily contain the highest numbers of methods or attributes. Neither is the nature of those classes suitable for use by subclasses. For example, String and Class are final and Hashtable and Vector are themselves subclasses of Dictionary and AbstractList, respectively.

#### 4.2.2 The Tomcat system

Table 5 shows, for the Tomcat system, the numbers of methods, fields and constructors and the depth in the inheritance hierarchy, for the top fifteen classes when ranked in descending order by their reverse aggregation PG values.

From Table 5, it is interesting that the mean inheritance depth of the fifteen classes is considerably smaller than in the JDK. Only one class, Node.Root extends any class other than java.lang.Object, compared with six in the JDK. It is also interesting to note that no class in the top fifteen when ranked by inheritance is present in either of the lists for top classes when ranked by normal or reverse aggregation PG values.

We would expect the depth values for both Tomcat and Ant to be smaller than those for the JDK, and this is borne out in the values from Table 5. A contributing factor to this is the role that the JDK class java.lang.Object plays - every class in the JDK extends class java.lang.Object by default.

The top fifteen classes were then extracted for normal aggregation PG values. Table 6 shows that six of these top fifteen could also be found in the fifteen reverse aggregation values from Table 5.

We note that the class JspC with 57 attributes (see Table 1) did not figure in either the highest normal or reverse aggregation PG values. Inspection of the raw data revealed it to be ranked 172 (reverse aggregation) and 73 (normal aggregation) from 321 classes analysed. Visual inspection revealed that most attributes of the class were Strings and ints. Also, the class with 168 methods (see Table 1), catalina.core.StandardContext, has 49 attributes and is ranked 37th by reverse aggregation PG, 14th by aggregation PG and 41st by inheritance. Although it does not fit in the top fifteen by reverse aggregation, it may be considered a candidate class for refactoring. The most obvious refactorings relate to a bad smell - namely Primitive Obsession .

With regards to conjecture C1, it is not true to say that classes with substantial amounts of aggregation (either normal or reverse), or indeed, with substantial amounts of methods and/or attributes tend to have the highest number of descendents (i.e., tend to be found near the root of an inheritance hierarchy).

With respect to refactoring of constructors, we note that only two classes of those listed in Table 5 have three or more constructors. This would imply less scope for refactoring by replacing multiple constructors with creation methods in this system. This constrasts with the JDK where there were six such classes.

#### 4.2.3 The Ant system

Table 7 shows values the the fifteen highest reverse aggregation values. None of these classes were in the top fifteen classes when sorted by aggregation PG values. Only one class, Task was present in the top fifteen when sorted by inheritance PG values. This class was ranked 30th in the aggregation list - just outside the top 1%. One possible explanation for the lack of overlap is that the inheritance hierarchy in Ant centers around the Task and ProjectComponent classes. This means that the other classes in Ant depend more upon inheritance and less upon aggregation and delegation.

One hypothesis to explain this behaviour is that aggregation is used as a surrogate for inheritance. The results from the JDK and Tomcat systems, where there was very little overlap between inheritance and aggregation PG values, and in particular from the Ant system, suggests that this may be the case. The hypothesis that delegation is used as a surrogate for multiple inheritance fits well with the observation that the distribution of interfaces follows a power-law . Whilst many classes implement a few interfaces, a few classes implement a large number of interfaces. Those that do, tend to delegate the responsibility for the methods of these interfaces to members of the same interface.

In the Ant system, we would consider the top ten classes found by ordering reverse aggregation PG value to be key classes. Of these, Project with 89 methods and 33 attributes stands out as Large Class and a possible candidate for refactoring. Hence, conjecture C1 would seem to be supported for the Ant system.

With regards to refactoring, only one class of those in Table 7 has three constructors. The rest are ineligible for refactoring by replacing multiple constructors with creation methods. Together with the result for Tomcat, this implies a significant difference between libraries and application systems; namely key classes in library systems tend to have a greater number of constructors.

## 5 Meta-analysis of refactorings

The focus of the research in this paper has been to identify candidate classes for refactoring - i.e. key classes. We have chosen specific properties to identify such classes. We now need to justify why, for example, we chose aggregation and inheritance relationships, together with the number of methods and attributes as those features (as opposed to any other class features).

To illustrate why, and as part of the research herein, a dependency diagram showing the relationships between the seventy-two refactorings outlined in Fowler et al.  was developed. From this meta-analysis emerged various core refactorings. By core refactorings, we mean those refactorings upon which a large number of other refactorings depend and are required in each of those refactorings. The three most important of these, in terms of the number of refactorings dependent on them, as established by the meta-analysis, are:

1. Extract Method \- which should be performed when the code body of a method is getting too long.
2. Move Field \- which should be performed when that field is being used by another class more than by the class in which it is defined.
3. Move Method \- which should be performed whenever a method is using features of another class more than those in which it is defined.

For both of the latter two core refactorings, and to a lesser extent extract method, aggregation plays a central role. If we wish to move a field, then the type of that field must be considered; it may be that of another class. If we want to move a method from one class to another, then all types of coupling in that method (including inheritance and aggregation) must be considered; removing the coupling from the source class may simplify it, but their loss needs to be replaced with appropriate code. Analysis of the mechanics of the two refactorings reveals inheritance to play a large role, in keeping with many other refactorings. For example, the Move Field refactoring requires as part of its mechanics that if the field is not declared as private, then all subclasses need to be checked for references to that field. The Move Method refactoring requires all subclasses and superclasses to be checked for references to that method.

We also believe that classes with large numbers of methods and attributes are likely to require application of these two refactorings at some point. Moreover, we refer to two common bad smells in classes identified in this paper from the key classes chosen - that of Large Class characterised by too many methods and Primitive Obsession characterized by overuse of primitive attributes and system classes. We therefore justify our choice of criteria for selection of key classes on this basis of these arguments.

## 6 Conclusions and Future Research

In this paper, we have described a technique by which key classes can be identified from three Java systems. A metric, potential gain, was adapted from use in Web search and navigation problems to evaluate the dependence which the system places on various classes. By computing values for potential gain on graphs for two forms of coupling, namely, inheritance and aggregation, we were able to identify candidate key classes. By combining these results with two other class features (number of class methods and class attributes) we were able to further isolate classes in possible need of refactoring.

Four principal results emerged from the research. Firstly, that metrics from other research domains can be adapted to aid developers and researchers in the refactoring process.

Secondly, there are substantial differences between each of the three systems investigated. The JDK system has key classes with more constructors. This is reflected in the mean number of constructors for classes within these systems - 1.242 for JDK against 1.073 for Ant and 1.069 for Tomcat.

Thirdly, there is also evidence to suggest that only the Ant system exhibits the expected properties of key classes according to conjecture C1. We therefore reject the hypothesis that key classes are consistently found at the base of inheritance hierarchies.

Finally, the analysis supports previous work on power-laws, which suggested that interfaces were commonly used as a surrogate for multiple implementation inheritance.

Future work will focus on two areas. Firstly, expanding the scope of what a key class is. The PG metric will be used to extract details relating to method parameters, method return types and interfaces. A second area of future research will investigate the potential for the key classes identified in this paper to be refactored. The research thus represents a first step in establishing the features of classes most eligible for refactoring.

## Appendix A Potential Gain

Formally, if the fraction of all possible paths in a directed graph, $`G=(N,E)`$, to a depth, $`d>0`$, which start from a node, $`nN`$ is given by

$$R_d(n)=\underset{yOut(n)}{}\frac{R_{d1}(y)}{_{jN}R_{d1}(j)}$$

where $`R_0=1`$ and $`N`$ is the set of nodes, $`E`$ is the set of edges and $`Out(i)=\{j|(i,j)E\}`$, then the Potential Gain of $`n`$ is given by

$$Pg(n)=\underset{k=1}{\overset{d_{max}}{}}R_k(x)f(x)$$

The constant, $`R_0=1`$, whilst not strictly accurate, is used to to ensure that $`\mathrm{log}Pg(n)>0`$ holds for all $`n`$. Two reasonable functions for $`f(d)`$ are the reciprocal function:

$$f(d)=reciprocal(d)=1/d$$

and the exponential decay function:

$$f(d)=decay(d)=\gamma ^d$$

where $`0<\gamma <1`$ is a constant.

The original justification for the use of these measures in the context of Web search was based on the assumption that the utility of browsing a page diminishes with the distance of the page from the starting URL. This assumption is consistent with experiments carried out on real Web data , and with studies showing that the probability of a user following a path of length $`n`$ decreases as $`n`$ increases .

The justification for these measures in the context of coupling is based on the fact that the influence that two classes have on each other diminishes as the distance (in terms of coupling) increases. This is a features of many patterns - such as Mediator and Facade  which reduce communication and dependencies by introducing “middle-men”.

The PG metric was collected automatically for each type of coupling for each class in the three systems, and the reciprocal function was used to compute the potential gain values.


text:
# Agterberg and Dodgson Reply

## Abstract

A reply to the Comment by Mineev and Champel.

In their Comment min02 , Mineev and Champel have argued that our results are incorrect because of the simultaneous neglect of both non-linear and non-local terms in our theory. We disagree with this statement and point out that while Ref. min02  argue that these terms should be included, they do not give a convincing physical argument that these terms are important at low fields. In the following, we explicitly show that non-linear and non-local terms can be safely ignored near $`H_{c1}`$ where our results are valid. Furthermore, our results agree with the assertion of Ref. min02  that there is no true $`A`$ to $`B`$ phase transition in the vortex phase. In Ref. agt02 , we refer to the diverging correlation length that describes the $`A`$ to $`B`$ transition in the Meissner phase.

To justify our claims, we give the missing steps between Eq. 2 and Eq. 3 of Ref. agt02 . The expression used for $`\mathrm{\Psi }_{}`$ in the $`A`$ phase was

$$(\xi _{A,}^2𝐃^2)\mathrm{\Psi }_{}(𝐫)=\frac{\stackrel{~}{\kappa }}{\kappa }(D_xD_y+D_yD_x)\mathrm{\Psi }_+(𝐫)$$
(1)
where $`D_i=_i+i2\pi A_i/\mathrm{\Phi }_0`$. Eq. 3 in Ref. min02  is found by setting the $`𝐃^2`$ operator to zero. This operator cuts off the divergence that Ref. min02  point out in their Eq. 3. We work with $`\stackrel{~}{\kappa }/\kappa 1`$ and within a London approach. Taking $`\mathrm{\Psi }_+(𝐫)=|\mathrm{\Psi }_+|e^{i\varphi }`$, $`\mathrm{\Psi }_{}=e^{i\varphi }\stackrel{~}{\mathrm{\Psi }}_{}`$, and defining the superfluid velocity as $`𝐯=\varphi +2\pi 𝐀/\mathrm{\Phi }_0`$ we get

$$(\xi _{}^2+𝐯^22i𝐯^2)\stackrel{~}{\mathrm{\Psi }}_{}(𝐫)=|\mathrm{\Psi }_+|\frac{\stackrel{~}{\kappa }}{\kappa }(i_xv_y+i_yv_x2v_xv_y).$$
(2)
This last equation shows that $`\stackrel{~}{\mathrm{\Psi }}_{}`$ will be proportional to $`v`$ (for small $`v`$). Consequently, for large vortex separations, the non-linear term in $`v`$ on the right hand side and the terms with $`v`$ on the left hand side of Eq. 2 can be ignored ($`ve^{r/\lambda }/\sqrt{r\lambda }`$ and the derivatives give a factor $`1/\lambda v`$ for large $`r`$). After Fourier transforming and using the Maxwell relation $`𝐯=\frac{2\pi \lambda ^2}{\mathrm{\Phi }_0}\times 𝐁`$, Eq. 2 becomes

$$\stackrel{~}{\mathrm{\Psi }}_{}(𝐪)=|\mathrm{\Psi }_+|\frac{\stackrel{~}{\kappa }}{\kappa }\frac{q_x^2q_y^2}{\xi _{A,}^2+𝐪^2}\frac{2\pi \lambda ^2}{\mathrm{\Phi }_0}B(𝐪).$$
(3)
Using this expression in the free energy and minimizing with respect to $`B(𝐪)`$ gives the novel London equation that forms the basis of our results

$$\left[1+\lambda ^2𝐪^2+\frac{\lambda ^2\stackrel{~}{\kappa }^2}{\kappa ^2}\frac{(q_x^2q_y^2)^2}{\xi _{A,}^2+𝐪^2}\right]B(𝐪)=0.$$
(4)
This gives a 4-fold symmetry to the structure of a flux line out to the distance $`\xi _{A,}`$, which diverges at the A$``$B transition temperature.

Using Eq. 3, $`|\mathrm{\Psi }_{}(𝐫)|\frac{\stackrel{~}{\kappa }}{\kappa }|\mathrm{\Psi }_+|\frac{2\pi \lambda ^2B(𝐫)}{\mathrm{\Phi }_0}`$ (take $`q_x^2q_y^2𝐪^2`$ in Eq. 3). This, with $`\beta _1|\psi _+|^2/\kappa =1/\xi _+^2`$ and $`2\pi \lambda ^2B(𝐫)/\mathrm{\Phi }_0=\lambda ^2/d^2`$, we find:

$$\frac{\beta _1|\mathrm{\Psi }_{}|^4}{\kappa |𝐃\mathrm{\Psi }_{}|^2}\frac{\stackrel{~}{\kappa }^2}{\kappa ^2}\frac{\lambda ^2}{\xi _{A,+}^2}\frac{\lambda ^4}{d^4}$$
(5)
where $`d`$ is the distance between vortices. Ref. min02  uses this result as their Eq. 5. $`H_{c1}`$ marks a second order phase transition between the Meisner and the vortex phase. At this transition $`B=0`$, which implies that as $`HH_{c1}`$, $`d\mathrm{}`$. Consequently, sufficiently near $`H_{c1}`$, the non-linear term can be safely neglected. Mineev and Champel argue that $`d`$ is cutoff when $`d\lambda /\sqrt{\mathrm{ln}(\lambda /\xi _{A,+})}`$. Presumably, this accounts for the very small range of $`H`$ over which $`B`$ goes from $`H_{c1}`$ to zero (or $`d`$ goes from $`\lambda `$ to $`\mathrm{}`$). Nevertheless, this range is experimentally accessible, and furthermore, for flat plate-like samples the demagnetization factors force the applied field $`H_{app}B`$, thus making the large $`d`$ limit even more accessible. It is also argued that when Eq. 6 in Ref.min02  is not satisfied, then non-local corrections must be included (these are terms that are $`O(𝐪^4)`$ or larger in Eq. 4). To counter this, we point out that Eq. 4 implies new physics at low fields (for small $`q`$ and large $`d`$) while non-local terms give new physics at high fields (for large $`q`$ and small $`d`$). In fact, we find non-local terms are not important when $`\stackrel{~}{\kappa }/\kappa >\frac{\xi _{A,+}}{\lambda }\frac{\lambda }{d}`$. Consequently, non-local terms also become negligible in the large $`d`$ limit.

DFA was supported by an award from Research Corporation. MJWD was supported by an EPSRC Advanced Fellowship AF/99/0725.


text:
# Neutralino annihilation 𝛾-rays from clumps and the LMC

## 1 Dark matter clumps

High resolution N-body simulations have revealed the survival of considerable substructure within galactic halos. Assuming these substructure clumps are composed of annihilating neutralinos, the flux $`F`$ of a clump at distance $`d`$ with a density distribution $`\rho (r)`$ is

$$F=\frac{1}{2}\frac{1}{4\pi d^2}\frac{N_\gamma \sigma v}{m_\chi ^2}_0^R\rho ^2(r)d^3r,$$
(1)
where $`\sigma v`$ is the thermally averaged annihilation cross section, $`m_\chi `$ is the neutralino mass, and $`N_\gamma `$ is the number of photons per annihilation with energy above an assumed energy threshold. For the density profile of the clumps we use the Navarro, Frenk, and White (1996) (NFW) and the Moore et al. (1998) profiles which were found to describe adequately the dark matter halos in simulations. As an upper limit to the degree of central concentration of the actual clump profile we also consider the singular isothermal sphere (SIS). The three profiles behave at the center as $`r^1`$, $`r^{1.5}`$, and $`r^2`$, respectively. The results for the profiles and for the minimum and maximum clump masses used, are shown in Fig. 1.

The SUSY models that give a 5-$`\sigma `$ (or more) detection using an Atmospheric Cherenkov Telescope (ACT) with effective area $`A_{eff}=10^8\mathrm{cm}^2`$, energy threshold $`E_{th}=50`$ GeV, and 100 hours of observation, are all the models that lie above the corresponding line for each density profile. The backgrounds used to calculate the noise are the hadronic and electronic cosmic ray shower contributions. Clearly, massive clumps appear to be easily detectable, regardless of profile; less massive clumps may be detectable, depending on their degree of central concentration. For more details see Tasitsiomi and Olinto (2002). There are some issues with respect to the ability of dark matter clumps to survive tidal disruption at distances small enough to yield easily detectable fluxes. In addition, there are some observational issues, given that the exact location of these clumps is not known. Thus, we focus on an object whose location is known, the LMC.

## 2 Flux from the LMC

We derive the density profile needed to calculate the $`\gamma `$-ray flux using rotation velocity data (Kim et al., 1998; Alves and Nelson, 2000). We fit the data using both the NFW ($`\rho _{NFW}`$), and the Hayashi et al. (2003) profile ($`\rho _H`$),

$$\rho _{NFW}=\frac{\rho _0}{r/r_s(1+r/r_s)^2},\rho _H=\frac{\rho _{NFW}}{1+(r/r_t)^3}.$$
(2)
Numerous observations indicate that the LMC is tidally stripped. The Hayashi et al. profile is a modification of the NFW which accounts for the tidal stripping that a halo may have undergone.

The fits are shown in the left panel of Fig. 2. In the right panel we present our results for the part of the SUSY parameter space that gives a 5-$`\sigma `$ detection for GLAST and a typical ACT; the instrument and observation parameters used to derive these results are also shown. Requiring that $`F_{LMC}F_{min}^{GLAST}=10^{10}`$ cm-2 s-1 we find that GLAST will be able to detect the signal for a significant part of the parameter space. This is true especially if the recently derived limit $`m_\chi <500`$ GeV is taken into account(Ellis et al., 2003). Assuming standard specifications, ACTs will not be able to probe any part of the parameter space (dotted line). Unless, large integration times (say, $`1`$ month) and effective areas (say, $`10^9`$ cm2) are used (dashed line). Note though that such integration times are fairly long for ACT observations, and that such large effective areas for an energy threshold $``$ 50 GeV are beyond the goals of existing and upcoming ACTs. These conclusions remain essentially the same for all the profiles used to model the LMC halo. The spectrum and its characteristic features, such as the cutoff at $`E=m_\chi `$ will be useful in identifying neutralino annihilation as the origin of the observed flux, especially in the case of GLAST where the prospects of detection are fairly good. The monochromatic lines produced by neutralino annihilation (e.g., the $`\gamma \gamma `$ line at $`E=m_\chi `$) would be excellent observational signatures if the cross sections for these processes were not suppressed(Tasitsiomi and Olinto, 2002).

EGRET has detected a flux of $`(14.4\pm 4.7)\times 10^8`$ photons $`(E>100\mathrm{MeV})`$ cm-2 $`s^1`$ from the LMC (Hartman et al., 1999). The emission due to neutralino annihilation can be anywhere from $`10^{13}`$ to $`10^9`$ photons $`(E>100\mathrm{MeV})\mathrm{cm}^2\mathrm{s}^1`$, depending on the neutralino parameters. The maximum possible flux is still $``$ 2 orders of magnitude less than the observed flux. This means that cosmic rays may be almost the exclusive source of the observed flux, as is often assumed. This is verified in the case of synchrotron emission from neutralino annihilation as well \[for the synchrotron and for more details on the $`\gamma `$-rays see Tasitsiomi, Gaskins, and Olinto (2003)\].

## 3 Conclusions

Dark matter clumps are in principle detectable, depending on the SUSY parameters, their distance and their degree of central concentration. The expected $`\gamma `$-ray flux from the LMC is, for most SUSY models, significantly smaller than what EGRET observed, verifying the usual assumption that cosmic rays are almost exclusively the origin of the detected flux. However, the flux is high enough to renter a large part of the SUSY parameter space accessible to GLAST; the detection of the signal by ACTs is highly unlikely.