text:
# S986 in M67: A Totally-Eclipsing Binary at the Cluster Turnoff

## 1 Introduction

M67 is probably the most thoroughly studied old open cluster in the Galaxy, thanks to its small distance from us. Typically quoted values for the cluster’s age ($`4\pm 0.5`$ Gyr; Dinescu et al. 1995) place it between the majority of known open clusters and the much older globular clusters. A determination of the mass of turnoff stars in this cluster would immediately provide us with a means of critically testing the validity of theoretical stellar evolution tracks, and checking models of various physical effects (such as convective overshooting and chemical diffusion). However, because M67 seems to have a high binary star fraction ($`50`$%; Fan et al. 1996) and there is definite evidence of strong dynamical interactions among cluster members (e.g. van den Berg et al. 2001; Sandquist et al. 2003), special care must be taken to ensure that the stars in such a binary have not been influenced by their environment.

Accurate photometry of the star S986 \[ID number from Sanders 1977; also known as F111 (Fagerholm 1906), and MMJ 5624 (Montgomery, Marschall, & Janes 1993)\] indicates that it falls at the turnoff in the cluster’s color-magnitude diagram (CMD), and proper motion studies identify it as a high probability ($`95\%`$) member of the cluster (Sanders, 1977; Girard et al., 1989; Zhao et al., 1993). S986 was previously identified as a single-lined spectroscopic binary by Mathieu, Latham, & Griffin (1990) with a period $`P=10.3386\pm 0.0006`$ d. Another interesting aspect of this system is the fact that the primary is on a circular orbit, indicating that tidal interactions between the stars have damped out any initial eccentricity. Over time, circularized main sequence binaries should be found at longer and longer periods as these interactions have sufficient time to work in wider and wider binaries (Mathieu & Mazeh 1988). S986 is the longest period circularized binary known in the cluster, and has probably required most of the cluster’s lifetime to become circularized.

As Mathieu et al. (1990) discuss, the secondary star in the binary must be at least about 2 magnitudes fainter than the primary at 5200 Å because it is a single-lined system. This indicates that the primary must contribute most of the light coming from the system, and therefore must be very close to the turnoff mass. A lower limit on the mass of the secondary is placed by the mass function for the binary: approximately 0.51 $`M_{\mathrm{}}`$in the case of S986 (assuming a primary mass of 1.25 $`M_{\mathrm{}}`$). Mathieu et al. stated that residuals of the orbital solution correlated with phase, and that a weak secondary cross correlation peak could be seen. Later Melo, Pasquini, & De Medeiros (2001) definitively detected a secondary peak in the cross-correlation function for the system, and measured a rotational velocity $`v\mathrm{sin}i=5.3\pm 0.6`$ km s-1. Careful examination of their cross correlation function indicates that the secondary peak is rather close to systematic radial velocity for M67.

As a byproduct of our project to monitor the partially-eclipsing blue straggler S1082 (period 1.0677978 d), we made extensive observations of the fields near the core of M67. The period of S986 has unfortunate aliases with the period of Earth’s rotation so that only one good apparition of the primary eclipse occurs during nighttime hours for a given observing site per month. This may explain the reason that eclipses have not previously been detected for this system. In the following sections, we describe our spectroscopic and photometric observations and reduction, the modeling of the light curve, and finally a discussion of the astrophysical meaning.

## 2 Spectroscopic Observations and Reduction

Because the spectroscopic observations of S986 provide very important constraints on the properties of the stars involved, we discuss them first. Our spectra were obtained at the Hobby-Eberly Telescope (HET) with the High Resolution Spectrograph (HRS) as part of normal queue scheduled observing for the period 2003-1. Our initial spectra were taken with the purpose of investigating the detection of a secondary peak in the cross-correlation function (CCF) by Melo et al. (2001). The spectra were taken with the three arcsecond science fiber with one sky fiber at 60,000 resolution. Exposure times were 540 seconds and yielded signal-to-noise (S/N) of 28 per resolution element. Follow-up spectra were taken with the two arcsecond science fiber with two sky fibers at 15,000 resolution. Exposure times were 1000 seconds and yielded S/N of 100 per resolution element. The spectra were reduced with IRAF111IRAF (Image Reduction and Analysis Facility) is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under contract with the National Science Foundation. ECHELLE scripts.

To determine radial velocities, we employed the IRAF task FXCOR to cross-correlate the reduced spectra against a solar spectrum taken with the same spectrograph. We confirmed the secondary peak in the CCF, and deblended the two components using Gaussian profiles with independent centroids and FWHM. The region used for determining the CCF was 5420 - 5720 Å for the 60,000 resolution spectra and 5960 - 6440 Å for the 15,000 resolution spectra. We find no discernible radial velocity variation of the secondary peak even though we observed the primary star at its maximum deviation from the cluster mean, and so we conclude that it is not the secondary component of the binary. Although we did examine the CCF for an even fainter third component, we did not find any significant traces. The secondary peak in the CCF sits at the cluster velocity (33.5 km s-1; Mathieu & Latham 1986), so that we conclude it is a cluster member: either a chance superposition or a third star on a wide orbit around the inner binary. We will denote the bright star in the binary as component Aa, and the third star as component B. The velocities measured for these two components are given in Table 1. The FWHM for the primary component Aa was found to be consistent with no detectable rotation ($`<`$ 10 km s-1) given our resolution, while the tertiary component B was found to have a small detectable rotational broadening ($`12\pm 4`$ km s-1).

The 60,000 resolution spectra were then shifted to the zero velocity for component Aa and combined to produce a higher S/N master spectrum that was used for an abundance analysis. The techniques employed here are nearly the same as those used in Shetrone & Sandquist (2000): equivalent widths were measured and abundances computed for each line. The temperature, surface gravity, and microturbulence were varied to minimize the slopes in the abundance vs. excitation potential and abundance vs. reduced equivalent width plots, and also to force ionization equilibrium. As an additional constraint, we limit the microturbulence to a value consistent with the formula given by Edvardsson et al. (1993) as a function of surface gravity and effective temperature to within their quoted errors (0.3 km s-1). One departure from the analysis of Shetrone & Sandquist (2000) was the use of a flux contribution variable that was divided into Aa’s EW to compensate for the contribution of the B spectral component (although the B component’s spectral lines were removed the flux contribution remains in our master Aa spectrum). The flux contribution variable was set by forcing the overall derived metallicity to be the same as the cluster mean ($`0.05`$; Shetrone & Sandquist 2000). Using this technique we found the primary component has a $`T_{eff}=6400\pm 50`$ K, $`\mathrm{log}g=4.25\pm 0.08`$, $`v_t=2.0\pm 0.1`$ km s-1and a flux contribution of $`0.885\pm 0.015`$. Because the “third light” contribution to the system flux by component B affects the observed depths of the eclipses described below, this measurement is an important ingredient in determining system parameters.

As part of the abundance analysis we detected a small Li line in the master spectrum of the Aa component having an equivalent width of 10 mÅ. This corresponds to $`\mathrm{log}N(Li)=2.11_{0.4}^{+0.2}`$ (see Fig.1), which is near the lower limit of 2.04 derived by Jones, Fischer, & Soderblom (1999) from data from García Lopéz, Rebolo, & Beckman (1988) (although without any corrections for the other components). In addition, our examination of O, s-process, and $`\alpha `$ elements for component Aa seem to be normal for main sequence or blue straggler stars (Shetrone & Sandquist 2000).

A master B spectrum was created by dividing the master Aa component into the individual spectra and then shifting the spectra to the rest velocity of the B component and combining. Unfortunately, because of the low S/N of the R=60,000 spectra, individual EW could not be measured from that master B spectrum, although the R=15,000 spectra had sufficient signal to identify the weak B component lines (typical EW = 8 mÅ ). An abundance analysis technique similiar to the one described above was applied to this B component but due to the small number of lines and high noise, we were forced to constrain the surface gravity using the assumption that the B component is a main sequence dwarf with a surface gravity and microturbulence appropriate to its $`T_{eff}`$. Using this technique we find that the B component has $`T_{eff}=5750\pm 200`$ K with a flux contribution of $`0.097\pm 0.025`$.

## 3 Photometric Observations and Reductions

All of the photometry for this study was taken at the 1 m telescope at the Mt. Laguna Observatory using a $`2048\times 2048`$ CCD on nights between December 2000 and November 2002. The nights of observations are given in Table 2. The photometry was in $`V`$ and $`I`$ bands with typical exposure times of 20 s (ranging between 15 and 60 s). Exposures were usually separated by about 2.5 minutes due to a relatively long readout time for the CCD.

Most of the details of the reduction are presented in other papers (Sandquist et al. 2003, Sandquist & Shetrone 2003), so we only briefly describe the reduction here. The object frames were reduced using overscan subtraction, bias frames, and flat fields. We conducted aperture photometry using the IRAF tasks DAOFIND and PHOT from the APPHOT package. In order to improve the accuracy of the relative photometry for the light curves, we used an ensemble photometry method similar to that described by Honeycutt (1992), iterating toward a consistent solution for photometric zeropoints for all frames and average magnitudes for all stars. Because frames typically had more than 300 measurable stars, the formal errors in the zero points ranged from around 0.003 to 0.007 mag, even with respect to night-to-night variations. Our $`V`$ and $`I`$ data (given relative to the median magnitude) are presented in Tables 3 and 4, respectively.

Our photometry has been calibrated using cluster frames taken during the same photometric night as a large number of Landolt (1992) and Stetson (2000) standard stars. We leave the details of the calibration and the comparison with previous datasets to a separate paper (Sandquist 2003).

## 4 Light Curves

Using our observations, we first updated the ephemeris of Mathieu et al. (1990). We have four sequences of photometric observations covering the ingress to primary eclipse, and one set of observations covering a portion of an egress. Because the system appears to be non-interacting, it is reasonable to expect that the period of the system has remained constant since the radial velocity observations tabulated by Mathieu et al. (1990). Times of eclipse contacts were determined by $`\chi ^2`$ minimization for the light curves relative to a trial model light curve. The times of observed first contact were HJD 2451959.83, 2452352.68, 2452383.69, and 2452662.83, while a fourth contact was observed at HJD 2452228.95, and one mid-eclipse point could be determined at HJD 2452600.96. This information was used in a simultaneous five-parameter fit to the combined dataset, where the parameters used were period $`P`$, epoch of primary eclipse $`T_0`$, system velocity $`\gamma `$, velocity semi-amplitude $`K_1`$, and phase half-width of primary eclipse $`\delta \varphi _p`$, assuming that the orbit has zero eccentricity. (Tests indicated that there was no measurable eccentricity.) The ephemeris for primary eclipse is

$$2445788.13(1)+10.33813(7)E$$

The numbers in parentheses indicate the uncertainty in the last digits of $`T_0`$ and $`P`$. The period is consistent (as expected) with the value of Mathieu et al. to within their quoted error. The velocity parameters also match the Mathieu et al. measurements very well since the same radial velocity data was used: $`\gamma =33.7\pm 0.1`$ km s-1and $`K_1=33.8\pm 0.2`$ km s-1. Finally, the fit constrains the width of primary eclipse to be $`\varphi _4\varphi _1=2\delta \varphi _p=0.0318\pm 0.0084`$. The reduced $`\chi ^2`$ for the overall fit was 1.32.

Phased data for the phases of primary and secondary eclipse are shown in Figs. 2, 3, and 4. The curvature of the light curve at second and third contacts (entering and exiting total eclipse) indicates that the primary eclipse is indeed due to a smaller star transiting across the face of a larger star.

### 4.1 Eclipse Depths

The eclipse depths provide robust information about the flux ratios of the two stars (although this is affected by the third star contributing to the system’s light). The overall depth of the primary eclipse is approximately 0.086 mag in $`V`$ and 0.078 mag in $`I`$, consistent with eclipse by a faint cool companion. The secondary eclipse was not detected in $`V`$, but we did find evidence of a decrease in brightness in $`I`$ on several nights. Data taken on the night with best atmospheric conditions (Jan. 17/18, 2003) provided our cleanest measurement of the depth of the secondary eclipse in $`I`$ ($`0.0111\pm 0.0011`$ mag) using the difference between the average magnitudes out of eclipse and in total eclipse. (It appears that there was a slight zero-point differences between the two nights, so a single reference level in $`I`$ was not used. This may have come about because the majority of our measurements in $`I`$ were taken during nights of eclipse.) We have a second measurement of the eclipse depth ($`0.0101\pm 0.0013`$ mag) from Apr. 20/21, 2003. Combining these measurements, we have a final value $`0.0107\pm 0.0008`$ mag.

### 4.2 Inclination and Radius Ratio

The ratio of the radii of the two stars is also a robust measurement because it depends on the timing of inner and outer contacts. Even with an uncertain contribution from third light, the radii of the two stars as determined from models are correlated, which means that the radius ratio is minimally affected. With fairly high time resolution observations, the times of contact can be measured directly. However, due to the duration of the eclipses, we have not yet been able to observe more than two eclipse contacts per night, so that some light-curve modeling is necessary. The phase width of the primary eclipse primarily constrains the inclination of the system and the sizes of the two stars relative to the orbital separation. Though we do not have a radial velocity curve for the secondary star, the resultant uncertainties in the mass ratio $`q`$ do not significantly affect the light curve fit. Because the secondary star acts like an almost perfectly black circular mask for most of the purposes of the modeling the primary eclipses, uncertainties introduced by error in the effective temperatures of the components are also minimal.

We have used the program NIGHTFALL222See http://www.lsw.uni-heidelberg.de/ rwichman/Nightfall.html for the program and a user manual (Wichmann 1998), which allows the use of model atmospheres (Hauschildt et al. 1999) and the modeling of physical effects such as detailed reflection. Detailed fitting of the primary eclipse requires some care in choosing the limb-darkening law. We have chosen to use a two-parameter square-root law based on literature comparisons with model atmospheres (Van Hamme, 1993; Díaz-Cordovés, Claret, & Giménez, 1995). Trials indicate that there are no observable indications of ellipsoidal shapes for the components, or significant reflection effects.

We computed light curve models covering a grid in the flux fraction of the primary $`f_{Aa}(V)`$, and the inclination of the binary $`I`$. We have imposed a constraint that the star contributing the third light should be on the single star fiducial line for M67. This assumption is not important for the determination of the inclination, but it has a slight effect on the radius determinations, especially if the star is itself an unresolved binary.

Although the uncertainty in the third light contribution has a significant effect on the uncertainties in the stellar radii, the uncertainty in the inclination dominates. The stellar radii are correlated with the inclination, with smaller inclinations requiring larger stellar radii to maintain the totality of the eclipse. Based on $`\chi ^2`$ minimization, the inclination is $`89\mathrm{°}_3^{+1}`$ where the lower limit is set by when the primary eclipse unavoidably becomes partial. A more stringent lower limit can be derived based on when differences between observations and the best-fit model become significant at the $`3\sigma `$ level. Based on this we quote an inclination $`i=89\mathrm{°}_2^{+1}`$. In Figs. 2 and 3, we compare our data with best fit models of two different inclinations. The model fits at lower inclination first begin to become unsatisfactory at the contact points ($`\varphi \pm 0.015`$ and $`\pm 0.01`$). Using this constraint, the ratio of the stellar radii is $`R_{Ab}/R_{Aa}=0.268_{0.004}^{+0.015}`$, where the uncertainties include the uncertainty in $`f_{Aa}(V)`$.

### 4.3 Photometry of the Primary Star and its Identity

The measured $`T_{eff}`$ and flux contribution from component B (providing the third light) are consistent with it being a normal main sequence star, based on comparison with Yi et al. (2001) and Girardi et al. (2000) isochrones using a distance modulus $`(mM)_V=9.72`$ (Sandquist 2003). In the following analysis, we will assume that component B is a normal main sequence star and a cluster member. Based on that, we can constrain the photometric properties of the primary star in the eclipsing binary (component Aa) without resorting to models. Our calibration of the ensemble photometry (Sandquist 2003) leads to a value for the sum of the three stars of $`V_{tot}=12.729`$ and $`(VI)_{tot}=0.643`$. Fig. 6 shows the $`VI`$ CMD for our deconvolution of component B from the total light of the eclipsing binary. The three linked points delineate the range of $`f_{Aa}(V)`$ allowed by our spectra. Our derived values are $`V_{bin}=12.86\pm 0.02`$ and $`(VI)_{bin}=0.619\pm 0.002`$.

Assuming again that the secondary component of the eclipsing binary Ab is a normal main sequence member of the cluster, we can deconvolve the (small) contribution of its light and get an estimate of the photometry of the primary alone. From the measured secondary eclipse depth and estimates of the relative flux contributions of components Aa and B, we find $`I_{Ab}I_{Aa}=4.85`$, or $`I_{Ab}=17.10_{0.08}^{+0.09}`$. Component Ab essentially does not contribute in $`V`$ (only about 0.004 mag), so that the corrected primary star values are $`V_{Aa}=12.86\pm 0.02`$ and $`(VI)_{Aa}=0.611\pm 0.003`$.

The deconvolved photometry values for the primary place it significantly to the blue of the most densely populated portions of the turnoff. We have 2152 observations of S986 in $`V`$ and 816 in $`I`$ and comparable numbers of observations of many other turnoff stars, so that error in the photometry relative to turnoff stars can be ruled out. If we are wrong in our assumption that component B is a normal main sequence star and it is actually an unresolved binary itself, this would make B redder and would require component Aa to be bluer still. The effective temperature measured from our spectra indicates that it is hotter than other stars that are clearly on the single-star sequence in the color-magnitude diagram.

We can compare our observations with theoretical models based on the spectroscopic $`T_{eff}`$ and $`\mathrm{log}g`$, thereby avoiding complications with color-effective temperature relations. In Fig. 5, we compare against the points that have the highest $`T_{eff}`$ on the theoretical isochrones of Girardi et al. (2000; hereafter Padova) and Yi et al. (2001; hereafter Y2). Component Aa is clearly inconsistent with the turnoff points from either set of isochrones, although neither set of isochrones provides an entirely satisfactory fit to the turnoff in the color-magnitude diagram. Although the chosen amount of convective overshooting in the stellar models does affect this, comparisons of the photometry with isochrones indicates that the amount of overshooting used in the models is slightly too large for M67 (Sandquist 2003). A smaller amount of overshooting would tend to increase the difference between component Aa and the model turnoffs. However, if component Aa is a blue straggler, then its surface gravity is likely to be larger than that of a single turnoff star (which would have evolved farther from the zero-age main sequence), consistent with the observations.

### 4.4 Secondary Star Properties

Although component Aa does not appear to be on the single-star sequence for M67, it is quite close, and we can use that fact to constrain the properties of the fainter component Ab. An absolute comparison with theoretical models for low-mass stars is impossible because the secondary has only been detected by its eclipse in $`I`$ band, data on the eclipsing binary does yet allow a complete solution of the system, and the calibrated color range for our photometry does not cover the color of component Ab. In addition, there are substantial model differences in the treatment of surface boundary conditions for low-mass models (Chabrier & Baraffe 1997), and theoretical uncertainties involving $`T_{eff}`$-color transformations (Houdashelt, Bell, & Sweigart 2000; Yi et al. 2001). In particular, there is strong evidence that models for low mass stars become too blue on the lower main sequence (von Hippel et al. 2002; Sandquist 2003).

The most direct and robust information about component Ab comes from the radius ratio $`R_{Ab}/R_{Aa}`$ and the eclipse depth in $`I`$ band $`\mathrm{\Delta }I=I_{Aa}I_{Aab}`$. The value of the eclipse depth, corrected for the third light contribution, is $`\mathrm{\Delta }I=0.0124\pm 0.0009`$. Thus, any model comparisons we make have to be relative to component Aa. We used isochrones of Girardi et al. (2000) and Yi et al. (2001) to find models that were consistent with the spectroscopic $`T_{eff}`$ and $`\mathrm{log}g`$ of component Aa in order to partially avoid uncertainties in $`T_{eff}`$-color transformations. These comparisons are shown in Fig. 7. The properties of component Aa are consistent with those of a star less than about 2.6 Gyr old, which can be understood if the star was formed in a collision of lower-mass stars earlier in the cluster’s history.

The isochrones allow us to get an estimate of the mass and radius of component Aa using the spectroscopic gravity, although details of the distribution of the star’s chemicals introduce some uncertainty (e.g. Sills et al. 1997). This is possible because of the mass-radius relation that exists along each isochrone, with younger isochrones giving the overall largest mass and radius values. We find nearly identical ranges from both the Yi et al. (2001) and Girardi et al. (2000) isochrones: $`M_{Aa}=1.24_{0.08}^{+0.19}`$ $`M_{\mathrm{}}`$ and $`R_{Aa}=1.37_{0.16}^{+0.24}`$ $`R_{\mathrm{}}`$. We can also estimate the $`I`$-band magnitude of component Aa using the isochrones and the spectroscopic measurements, and again we find similar ranges from the two sets of isochrones ($`M_{I,Aa}=3.08_{0.28}^{+0.20}`$). The primary source of uncertainty in all three estimates is from the measurement of the surface gravity.

In order to make the final comparison with our observed values from the eclipse analysis, we use the models of Baraffe et al. (1998; hereafter BCAH) since they currently provide the best fit to the lower main sequence of M67 (although they begin to become too blue by $`V17`$; Sandquist 2003). The BCAH models are also the ones that have been most thoroughly compared to observed colors of low-mass stars, and seem to agree well down to $`T_{eff}3600`$ K (BCAH). The isochrones of the Y2 and Padova groups employ different tables to convert from $`T_{eff}`$ to color, and get widely different results. The BCAH models most nearly agree with those of the Y2 or Padova groups near solar mass, where the agreement is forced by calibration to the Sun. We made a correction of 0.077 mag in $`I`$ to account for zero-point differences between the Padova and Y2 isochrones and those of BCAH for solar-mass models. Such a correction is justified by current uncertainties in the zero-points of $`T_{eff}`$-color relations for different sets of isochrones. Our final comparison is shown in Fig. 8. The figure indicates that component Ab is smaller and/or brighter in $`I`$ band than the models predict. The uncertainties in the spectroscopic properties of component Aa mean, however, that the inconsistency is only at a little over $`1\sigma `$.

A reasonable estimation of the effective temperature of star Ab is next to impossible, partly because there are not enough well-measured cool dwarf stars to set an empirical $`T_{eff}`$-color relation or calibrate a theoretical one. In addition, there is some evidence that the giant stars typically used to calibrate relations at low temperatures may follow a systematically different relation than the dwarfs (Houdashelt et al., 2000). The depth of secondary eclipse is sensitive to $`T_{eff}`$, but the interpretation requires the use of updated stellar atmosphere models because of strongly non-grey nature of the atmospheres of low-mass stars (e.g. Chabrier & Baraffe 1997). From the BCAH models and a distance modulus $`(mM)_I=9.67\pm 0.03`$ (Sandquist 2003), we roughly estimate that $`T_{eff}=3700\pm 200`$ K, but we again note that this is near where their models start to diverge from color-magnitude diagrams of field dwarfs. This corresponds to a mass $`M_{Ab}0.52`$ $`M_{\mathrm{}}`$, so that $`q0.41`$.

## 5 Discussion and Conclusions

We have presented the first evidence of total eclipses (including eclipses of the faint secondary star in $`I`$ band) in the system S986 using an extensive series of photometric observations. We have taken high-resolution spectroscopy of the system and have identified a third star that contributes to the light of the system. The third star appears to be a cluster member, but may or may not be physically associated with the eclipsing binary. The results of our analysis are given in Table 5.

The detailed analysis of our spectroscopy and photometry for the S986 indicates that the primary star in the eclipsing binary (component Aa) is a star that is slightly (but significantly) hotter than the turnoff of the cluster. Two stellar explanations for this exist. One possibility is that the primary is a normal main sequence star that is in a relatively short lived phase of its evolution. The gap in the CMD of M67 with $`12.85<V<13.1`$ corresponds to a rapid phase before and after core hydrogen exhaustion during which the central convection zone disappears and a shell fusion source is established. The size of the gap depends on exactly when the evolutionary timescale is small enough that few or no stars are likely to be found in the phase, given the total population of stars in the cluster. Several other members of M67 exist in the same portion of the CMD as component Aa (S489, S602, S610, S615, S1271, S1503, and S1575; Sandquist 2003), which may lend some credence to this idea. Because turnoff stars in clusters of M67’s age have small but significant convective cores, they put difficult constraints on the theory of convective overshooting that have not been satisfied as yet. Once they are, the CMD position of component Aa can be re-evaluated.

A second explanation is that component Aa is a blue straggler. It is somewhat difficult to piece together a scenario that can explain the system’s orbital and photometric properties, but the light curve analysis makes it clear that component Ab is a relatively normal main sequence star and that the two stars are completely detached. Component B is also consistent with being a normal main sequence star. As such, scenarios involving any kind of mass transfer to star Aa are unlikely.

Our detection of Li in component Aa is consistent with the abundances of stars on the edge of the “Li gap” (centered at $`T_{eff}6700`$ K). On the other hand, there has not been a detection of Li in a blue straggler in M67 to date. If component Aa was created in the merger of two stars, the detection of Li would require that the more massive star to have retained a substantial amount of surface lithium in the time before the collision, that little mass was lost from the surface of the star in the merger (see Lombardi et al. 2002 for a discussion), and the remaining Li was not substantially depleted after the merger. For a star with $`T_{eff}=6400`$ K, the Li depletion timescale seems to be several Gyr (Jones et al. 1999), comparable to the age of M67, so post-merger depletion is probably not large enough to matter here. For a star to have retained surface Li for approximately 1.5 Gyr, $`T_{eff}`$ could be as low as about 5700 K based on observations of NGC 752 (Hobbs & Pilachowski 1986). This corresponds to a mass of approximately 1 $`M_{\mathrm{}}`$, meaning that the other star involved would had to have been low in mass ($`0.25`$ $`M_{\mathrm{}}`$). A blue straggler scenario is possible with the detection of Li, but the possible formation routes are fairly tightly constrained.

One of the most interesting details about the S986 system is that it is circularized in spite of its relatively long period and the main sequence nature of its component stars. The tidal circularization places a lower limit on the age of the binary, although the limit is not likely to be very strong given the state of our knowledge of tidal interactions (e.g. Terquem et al. 1998) and that S986 is often used to constrain the strength of tidal dissipation (after assuming an age for M67). Our best blue straggler scenario for the S986 system involves a collision of relatively low-mass stars during a close single-binary or binary-binary interaction. During the the interaction two of the stars could merge to form the primary star, while a third star could be captured onto a relatively close orbit. Because strong multiple star interactions typically produce eccentric orbits, a period of circularization would still be necessary, which still should have taken most of the age of the cluster. However, this would also have the benefit of explaining why component Aa is so close to the cluster turnoff: the merger would occur before either of the stars had done much nuclear processing of their material. (If the input stars were both of low mass, this would also help.) Our isochrone comparisons imply that the merger probably happened more than 1.4 Gyr ago. Component B may or may not be associated with the creation of the system in this scenario. Long-term radial velocity monitoring would help answer that question.

Whether or not component Aa is a normal main sequence star or a blue straggler, further study of the binary is worthwhile because it provides a means of simultaneously comparing the properties of a low-mass main sequence star and a star with slightly more mass than the Sun. Because the structures of these two types of stars are rather different, theoretical models often require assumptions that can result in systematic errors when isochrones are computed. In the case of an open cluster like M67, the degree of convective core overshooting is an excellent example since it substantially affects the evolution of stars slightly more massive than the Sun (having relatively small but significant convective cores), but has negligible effect on lower main sequence models, which have surface convection zones of large extent. The surface boundary condition is also a potential source of problems because it sets the interior structure of fully- or mostly-convective stars while having a minor effect on models of stars with mass near that of the Sun (Chabrier & Baraffe 1997). M67 provides an excellent testing ground for these kinds of stellar physics issues.

E.L.S. would like to thank the director of Mount Laguna Observatory (P. Etzel) for generous allocations of telescope time to constrain the ephemeris of this binary system, L. Detweiler for allowing us to use one of his nights to observe an eclipse, and P. Etzel and D. Latham for useful conversations. This research was supported in part by the National Science Foundation under Grant No. AST-0098696 to E.L.S.


format: text

text:
# 1 The T-matrix method

## 1 The T-matrix method

The T-matrix method in wave scattering involves writing the relationship between the wave incident upon a scatterer, expanded in terms of orthogonal eigenfunctions,

$$U_{\mathrm{inc}}=\underset{n}{\overset{\mathrm{}}{}}a_n\psi _n^{(\mathrm{inc})},$$
(1)
where $`a_n`$ are the expansion coefficients for the incident wave, and the scattered wave, also expanded in terms of orthogonal eigenfunctions,

$$U_{\mathrm{scat}}=\underset{k}{\overset{\mathrm{}}{}}p_k\psi _k^{(\mathrm{scat})},$$
(2)
where $`p_k`$ are the expansion coefficients for the scattered wave, is written as a simple matrix equation

$$p_k=\underset{n}{\overset{\mathrm{}}{}}T_{kn}a_n$$
(3)
or, in more concise notation,

$$𝐏=\mathrm{𝐓𝐀}$$
(4)
where $`T_{kn}`$ are the elements of the T-matrix. The T-matrix method can be used for scalar waves or vector waves in a variety of geometries, with the only restrictions being that the geometry of the problem permits expansion of the waves as discrete series in terms of orthogonal eigenfunctions, that the response of the scatterer to the incident wave is linear, and that the expansion series for the waves can be truncated at a finite number of terms. In general, one calculates the T-matrix, although it is conceivable that it might be measured experimentally.

The T-matrix depends only on the particle—its composition, size, shape, and orientation—and is independent of the incident field. This means that for any particular particle, the T-matrix only needs to be calculated once, and can then be used for repeated calculations. This is a significant advantage over many other methods of calculating scattering where the entire calculation needs to be repeated . Some cases provide even more efficiency: if the waves are expanded in spherical functions, the averaging of scattering over various orientations of the particle compared to the direction of the incident wave can be performed analytically .

In the spherical geometry of elastic light scattering by a particle contained entirely within some radius $`r_0`$, the eigenfunction expansions of the fields are made in terms of vector spherical wavefunctions (VSWFs) :

$`𝐌_{nm}^{(1,2)}(k\mathrm{r})`$ $`=`$ $`N_nh_n^{(1,2)}(kr)𝐂_{nm}(\theta ,\varphi )`$ (5)
$`𝐍_{nm}^{(1,2)}(k\mathrm{r})`$ $`=`$ $`{\displaystyle \frac{h_n^{(1,2)}(kr)}{krN_n}}𝐏_{nm}(\theta ,\varphi )+`$ (6)
$`N_n\left(h_{n1}^{(1,2)}(kr){\displaystyle \frac{nh_n^{(1,2)}(kr)}{kr}}\right)𝐁_{nm}(\theta ,\varphi )`$
where $`h_n^{(1,2)}(kr)`$ are spherical Hankel functions of the first and second kind, $`N_n=1/\sqrt{n(n+1)}`$ are normalisation constants, and $`𝐁_{nm}(\theta ,\varphi )=𝐫Y_n^m(\theta ,\varphi )`$, $`𝐂_{nm}(\theta ,\varphi )=\times \left(𝐫Y_n^m(\theta ,\varphi )\right)`$, and $`𝐏_{nm}(\theta ,\varphi )=\widehat{𝐫}Y_n^m(\theta ,\varphi )`$ are the vector spherical harmonics , and $`Y_n^m(\theta ,\varphi )`$ are normalised scalar spherical harmonics. The usual polar spherical coordinates are used, where $`\theta `$ is the co-latitude measured from the $`+z`$ axis, and $`\varphi `$ is the azimuth, measured from the $`+x`$ axis towards the $`+y`$ axis.

$`𝐌_{nm}^{(1)}`$ and $`𝐍_{nm}^{(1)}`$ are outward-propagating TE and TM multipole fields, while $`𝐌_{nm}^{(2)}`$ and $`𝐍_{nm}^{(2)}`$ are the corresponding inward-propagating multipole fields. Since these wavefunctions are purely incoming and purely outgoing, each has a singularity at the origin. Since fields that are free of singularities are of interest, it is useful to define the singularity-free regular vector spherical wavefunctions:

$`\mathrm{𝐑𝐠𝐌}_{nm}(k\mathrm{r})`$ $`=`$ $`\frac{1}{2}[𝐌_{nm}^{(1)}(k\mathrm{r})+𝐌_{nm}^{(2)}(k\mathrm{r})],`$ (7)
$`\mathrm{𝐑𝐠𝐍}_{nm}(k\mathrm{r})`$ $`=`$ $`\frac{1}{2}[𝐍_{nm}^{(1)}(k\mathrm{r})+𝐍_{nm}^{(2)}(k\mathrm{r})].`$ (8)
Since the spherical Bessel functions $`j_n(kr)=\frac{1}{2}(h_n^{(1)}(kr)+h_n^{(2)}(kr))`$, the regular VSWFs are identical to the incoming and outgoing VSWFs except for the replacement of the spherical Hankel functions by spherical Bessel functions.

Since the incident field, in the absence of a scatterer, is singularity-free, the expansion

$$𝐄_{\mathrm{inc}}(\mathrm{r})=\underset{n=1}{\overset{\mathrm{}}{}}\underset{m=n}{\overset{n}{}}a_{nm}^{(3)}\mathrm{𝐑𝐠𝐌}_{nm}(k\mathrm{r})+b_{nm}^{(3)}\mathrm{𝐑𝐠𝐍}_{nm}(k\mathrm{r})$$
(9)
is generally used for the incident field. Alternatively, the purely incoming part of the incident field can be used:

$$𝐄_{\mathrm{inc}}(\mathrm{r})=\underset{n=1}{\overset{\mathrm{}}{}}\underset{m=n}{\overset{n}{}}a_{nm}^{(2)}𝐌_{nm}^{(2)}(k\mathrm{r})+b_{nm}^{(2)}𝐍_{nm}^{(2)}(k\mathrm{r}).$$
(10)
In both cases, the scattered field is

$$𝐄_{\mathrm{scat}}(\mathrm{r})=\underset{n=1}{\overset{\mathrm{}}{}}\underset{m=n}{\overset{n}{}}p_{nm}^{(1)}𝐌_{nm}^{(1)}(k\mathrm{r})+q_{nm}^{(1)}𝐍_{nm}^{(1)}(k\mathrm{r}).$$
(11)
The two sets of expansion coefficients for the incident/incoming field are related, since $`a_{nm}^{(3)}=2a_{nm}^{(2)}`$ and $`b_{nm}^{(3)}=2b_{nm}^{(2)}`$. However, the scattered/outgoing field expansion coefficients will differ, as will the T-matrix. Using the regular expansion, the T-matrix in the absence of a scatterer is a zero matrix, while using the incoming field expansion, the no-scatterer T-matrix is the identity matrix. The two expansions are essentially the same—the only difference is that the incident wave in the incident/scattered wave expansion includes part of the outgoing wave. T-matrices for the two expansions only differ by the identity matrix, so $`𝐓^{(\mathrm{in}/\mathrm{out})}=2𝐓^{(\mathrm{inc}/\mathrm{scat})}+𝐈`$. The incident/scattered formulation is much more commonly used; the incoming/outgoing formulation gives simpler results for the transport of momentum and angular momentum (that is, optical force and torque) by the field. It should be note that for plane wave illumination, for which the VSWF expansion is non-terminating, the incident/scattered formulation gives a scattered wave expansion that converges over all space, while the incoming/outgoing expansion, strictly used, would give an non-terminating, non-convergent outgoing field expansion. For focussed beam illumination with a finite VSWF expansion, the incoming/outgoing expansion directly gives the total outgoing field that would be experimentally measured. Since conversion from one formulation to the other is simple, either can be readily used for calculation of fields, forces, scattering matrices, or for orientation averaging.

In practice, the field expansions and the T-matrix are terminated at some $`n=N_{\mathrm{max}}`$. For the case of a scatterer that is contained within a radius $`r_0`$, $`N_{\mathrm{max}}kr_0`$ is usually adequate, but $`N_{\mathrm{max}}=kr_0+3\sqrt[3]{kr_0}`$ is advisable if higher accuracy is needed . Although we assume in this paper (as is usually the case) that the incident and scattered wave expansions are terminated at the same $`N_{\mathrm{max}}`$ (giving a square T-matrix), this is not necessary. It should be noted that convergence of the expansion of the incident field is not a necessary condition for the T-matrix method to be useful—indeed, for the most common application, scattering of an incident plane wave, the incident field expansion does not converge over all space. However, it does converge within the radius $`r_0`$—which is the part of the field that can affect the scattering particle—and therefore, the field expansions and the T-matrix can be truncated at a finite $`N_{\mathrm{max}}`$.

For the case of plane wave scattering, the plane wave expansion formula is useful:

$$a_{nm}=4\pi \mathrm{i}^nN_n𝐂_{nm}^{}𝐄_0,b_{nm}=4\pi \mathrm{i}^{n1}N_n𝐁_{nm}^{}𝐄_0.$$
(12)
The main case of interest for non-plane wave incident illumination is that of focussed beams. A variety of methods can be used, such as plane wave expansion , the localised approximation , or the point-matching method .

The only remaining requirement is that the T-matrix be calculated. This requires essentially a complete calculation of the scattering properties of the scatterer. This is almost universally done using the extended boundary condition method (EBCM), originally developed by Waterman , which is so strongly linked with the T-matrix method that the terms “EBCM” and “T-matrix method” are often used interchangeably. In the next section, we consider some general principles involved in the calculation of the T-matrix, and show that an alternative method—column-by-column calculation using the point-matching method (PMM)—is computationally feasible and simply implemented for homogeneous isotropic particles devoid of symmetry.

Lastly, before we continue to consider calculation of T-matrices in more detail, we can note that while the incident and scattered fields are usually expanded in terms of VSWFs, other sets of eigenfunctions, such as cylindrical wavefunctions (for scatterers of infinite length in one dimension), or a Floquet expansion (planar periodic scatterers), are more appropriate for other geometries. There is no requirement that the modes into which the incident and scattered fields are expanded be the same, or even similar. In all of these cases, the T-matrix method remains applicable.

## 2 Calculating the T-matrix

If the field expansions and T-matrix are truncated at some $`N_{\mathrm{max}}`$, there are $`N_T=2N_{\mathrm{max}}(N_{\mathrm{max}}+2)`$ expansion coefficients for each of the incident and scattered fields, and the T-matrix is $`N_T\times N_T`$. Since $`N_{\mathrm{max}}`$ is proportional to the radius enclosing the particle, $`r_0`$, the number of expansion coefficients is proportional to $`r_0^2`$, and the number of elements in the T-matrix is proportional to $`r_0^4`$. This can be used to obtain an estimate of the scaling of computational time for different methods of calculation.

### 2.1 The extended boundary condition method

In principle, any method of calculating scattering by the particle can be used to calculate the T-matrix. However, the method of choice is almost universally the EBCM . In the EBCM, the internal field within the particle is expanded in terms of regular VSWFs. Therefore, the method is restricted to homogeneous and isotropic particles. Rather than considering the coupling of the incident and scattered fields directly, the coupling between the incident and internal (the $`\mathrm{𝐑𝐠𝐐}`$ matrix), and scattered and internal fields (the $`𝐐`$ matrix) is calculated, and the T-matrix found from these ($`𝐓=\mathrm{𝐑𝐠𝐐𝐐}^1`$). The $`\mathrm{𝐑𝐠𝐐}`$ and $`𝐐`$ matrices are the same size as the T-matrix, with $`O(N_{\mathrm{max}}^4)`$ elements. The elements of these matrices are found by integrating over the surface of the scatterer, an operation requiring $`O(N_{\mathrm{max}}^2)`$ time per element, so the calculation of the $`\mathrm{𝐑𝐠𝐐}`$ and $`𝐐`$ matrices is expected to require $`O(N_{\mathrm{max}}^6)`$ computational time. The actual calculation of the T-matrix, if direct inversion is naïvely used, takes $`O(N_{\mathrm{max}}^6)`$ time. In practice, the calculation of the $`\mathrm{𝐑𝐠𝐐}`$ and $`𝐐`$ matrices dominates the computational time .

From this, it can be seen that the EBCM can be expected to be very slow for large particles. However, most applications of the EBCM have been for the special case of scattering particles rotationally symmetric about the $`z`$ axis. In this case, scattered modes $`𝐌_{n^{}m^{}}^{(1)}`$ and $`𝐍_{n^{}m^{}}^{(1)}`$ only couple to incident modes $`\mathrm{𝐑𝐠𝐌}_{nm}`$ and $`\mathrm{𝐑𝐠𝐍}_{nm}`$ if $`m^{}=m`$, greatly reducing the number of matrix elements that need to be calculated, and the surface integral over the particle surface reduces to a one-dimensional integral over $`\theta `$, since the azimuthal integration over $`\varphi `$ can be simply done analytically . This results in a great improvement in performance, and, in terms of computational time, EBCM is clearly the method of choice for axisymmetric particles. Numerical problems do occur when the scatterer is highly non-spherical. The discrete sources method is designed to overcome these problems . For the even more symmetric case of a spherical scatterer, the scattered and incident modes only couple if $`n^{}=n`$ and $`m^{}=m`$, the T-matrix becomes diagonal, and all of the integrals can be performed analytically, and Mie’s solution to scattering by a sphere  is simply obtained.

In a similar manner, scatterers with point-group rotational symmetry allow significant improvement of the computational time required through exploitation of the symmetry .

Methods have also been developed to calculate T-matrices for clusters of particles and for layered particles .

The efficiency of the EBCM for the calculation of the T-matrix is such that alternative methods need only be considered if the EBCM is inapplicable (such as when the particle in inhomogeneous or anisotropic), numerical difficulties are encountered using the EBCM (such as for extremely non-spherical particles), or if the scattering particle has no symmetries that can be used to optimise the computation of the T-matrix.

### 2.2 Methods other than the EBCM

Methods other than the EBCM can be used to calculate the T-matrix. In general, one would calculate the scattered field, given a particular incident field. The most direct way in which to use this to produce a T-matrix is to solve the scattering problem when the incident field is equal to a single spherical mode—that is, a single VSWF such as $`𝐄_{\mathrm{inc}}(\mathrm{r})=\mathrm{𝐑𝐠𝐌}_{11}(k\mathrm{r})`$, $`𝐄_{\mathrm{inc}}(\mathrm{r})=\mathrm{𝐑𝐠𝐍}_{11}(k\mathrm{r})`$, $`𝐄_{\mathrm{inc}}(\mathrm{r})=\mathrm{𝐑𝐠𝐌}_{21}(k\mathrm{r})`$, etc, and repeat this for all VSWFs that need to be considered (up to $`n=N_{\mathrm{max}}`$). The expansion coefficients for the scattered field can be found in each case, if necessary, by using the orthogonal eigenfunction transform (the generalised Fourier transform), and each scattering calculation gives a single column of the T-matrix.

Therefore, the calculation of a T-matrix requires that $`2N_{\mathrm{max}}(N_{\mathrm{max}}+2)`$ separate scattering problems are solved. The provides a criterion for deciding whether it is desirable to calculate a T-matrix: if more than $`2N_{\mathrm{max}}(N_{\mathrm{max}}+2)`$ scattering calculations will be performed, then it is more efficient to calculate the T-matrix and use this for the repeated calculations than it is to use the original scattering method repeatedly. Repeated calculations are expected if orientation averaging is to be carried out, or if inhomogeneous illumination is to be considered, such as, for example, scattering by focussed beams, where there are generally 6 degrees of freedom, namely the three-dimensional position of the scatterer within the beam, and the three-dimensional orientation of the scatterer. Even if only a modest number of points are considered along each degree of freedom, the total number of scattering calculations required rapidly becomes very large, and even if the T-matrix takes many hours to calculate, the total time saved by doing so can make an otherwise computationally infeasible problem tractable.

Volume methods are of interest, since they can readily be used for inhomogeneous or anisotropic particles. The two most likely candidates are the finite-difference time-domain method (FDTD)  and the discrete dipole approximation (DDA). In FDTD, the Maxwell equations are discretised in space and time, and, beginning from a known initial state, the electric and magnetic fields at each spatial grid point are calculated for successive steps in time. The number of grid points required is $`O(N_{\mathrm{max}}^3)`$ for three-dimensional scattering, and $`O(N_{\mathrm{max}})`$ time steps required, so FDTD solutions scale as $`O(N_{\mathrm{max}}^4)`$. Therefore, calculation of the T-matrix using FDTD should scale as $`O(N_{\mathrm{max}}^6)`$, which is the same scaling as the EBCM. However, the grid required must be closely spaced compared to the wavelength, and the space outside the scatterer must also be discretised, making FDTD substantially slower than EBCM, especially for smaller particles. However, FDTD is an extremely general technique, and has potential as a method for the calculation of T-matrices.

We should add that there is an additional consideration that makes FDTD potentially attractive as a method for calculating the T-matrix: FDTD does not assume that the incident wave is monochromatic. Consider the case when the illumination is a brief pulse with a Gaussian envelope. The frequency spectrum of the incident wave is Gaussian, and the scattering of a range of frequencies can be found by taking the Fourier transform of the scattered field . Even if we are not interested in other than monochromatic illumination, we will frequently be interested in scattering by size distributions of particles. Since varying the frequency for a particular particle is equivalent to varying the size of the particle for a fixed incident frequency, the T-matrices for a range of particle sizes can be calculated simultaneously.

The other major volume method for computational scattering, the discrete dipole approximation (DDA), also known as the coupled-dipole method, has been recently applied to the calculation of the T-matrix by Mackowski , who obtained good results, with reasonable computational efficiency using a moment method to solve the DDA system of equations. DDA lacks the main disadvantages of FDTD, namely the need to discretise space outside the particle, and the need to implement suitable boundary conditions to prevent non-physical reflections from the boundary of the computational domain. Mackowski’s method scales as $`O(N_{\mathrm{max}}^7)`$ for large $`N_{\mathrm{max}}`$. There is no need to discuss his method in detail here, and the interested reader is referred to his recent description of the method .

Finally, we consider the point-matching method. Like the T-matrix method and the EBCM, the point-matching method involves expansion of fields in terms of VSWFs. In the point-matching method, the internal field within the scatterer and the scattered field are expanded as series of VSWFs, and the incident, internal, and scattered fields are matched at points on the particle surface, using the usual boundary condition of continuity of tangential components of the electric and magnetic fields. This gives a system of equations from which the unknown expansion coefficients of the internal and scattered fields can be found. Typically, enough points are used for matching the fields so as to give an overdetermined system of equations, which is then solved in a least-squares sense. Solving the system of $`O(N_{\mathrm{max}}^2)`$ unknowns using direct matrix inversion can be expected to be an $`O(N_{\mathrm{max}}^6)`$ problem, with the result that the total computational time is $`O(N_{\mathrm{max}}^8)`$. In practice, faster methods can be used, and our results indicate a performance of about $`O(N_{\mathrm{max}}^7)`$ for our implementation.

The point-matching method is an attractive candidate since a T-matrix implementation will generally include routines to calculate VSWFs, making the implementation of a point-matching T-matrix calculator simple. The only further requirement is a routine for solving overdetermined linear systems. Since the scattered field is calculated in terms of VSWFs, the conversion of the results of a single PMM calculation to a T-matrix column is trivial.

Naturally, multiple expansion methods (the generalised multipole technique, or the multiple multipole method) can be used. Since multipole methods exist for anisotropic media , the method can be used for anisotropic scatterers.

Our implementation of the point-matching method, and its performance, is discussed in the next section.

## 3 Point-matching method

Our implementation of the PMM T-matrix calculation uses an incoming/outgoing field expansion (equations (10) and (11)), rather than the usual incident/scattered wave expansion (equations (9) and (11)), and the internal field is expanded in terms of regular VSWFs:

$`𝐄_{\mathrm{inc}}(\mathrm{r})`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{N_{\mathrm{max}}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}a_{nm}𝐌_{nm}^{(2)}(k\mathrm{r})+b_{nm}𝐍_{nm}^{(2)}(k\mathrm{r}),`$ (13)
$`𝐄_{\mathrm{scat}}(\mathrm{r})`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{N_{\mathrm{max}}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}p_{nm}𝐌_{nm}^{(1)}(k\mathrm{r})+q_{nm}𝐍_{nm}^{(1)}(k\mathrm{r}),`$ (14)
$`𝐄_{\mathrm{int}}(\mathrm{r})`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{N_{\mathrm{max}}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}c_{nm}\mathrm{𝐑𝐠𝐌}_{nm}(k\mathrm{r})+d_{nm}\mathrm{𝐑𝐠𝐍}_{nm}(k\mathrm{r}).`$ (15)
We use this particular expansion since we are interested in calculating optical forces and torques within optical traps  and this results in simpler expressions for these quantities.

We considered a single scatterer, centred on the origin, contained entirely within a radius $`r_0`$, and with a surface specified by a function of angle:

$$r=r(\theta ,\varphi )$$
(16)
The boundary conditions—matching the tangential fields on the surface of the scatterer—are

$`\widehat{𝐧}\times (𝐄_{\mathrm{inc}}(\mathrm{r})+𝐄_{\mathrm{scat}}(\mathrm{r}))`$ $`=`$ $`\widehat{𝐧}\times 𝐄_{\mathrm{int}}(\mathrm{r}),`$ (17)
$`\widehat{𝐧}\times (𝐇_{\mathrm{inc}}(\mathrm{r})+𝐇_{\mathrm{scat}}(\mathrm{r}))`$ $`=`$ $`\widehat{𝐧}\times 𝐇_{\mathrm{int}}(\mathrm{r}),`$ (18)
where $`\widehat{𝐧}`$ is a unit vector normal to the surface of the particle.

The magnetic fields are given by expansions similar to those for the electric fields:

$`𝐇_{\mathrm{inc}}(\mathrm{r})`$ $`=`$ $`{\displaystyle \frac{1}{k_{\mathrm{medium}}}}{\displaystyle \underset{n=1}{\overset{N_{\mathrm{max}}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}a_{nm}𝐍_{nm}^{(2)}(k\mathrm{r})+b_{nm}𝐌_{nm}^{(2)}(k\mathrm{r}),`$ (19)
$`𝐇_{\mathrm{scat}}(\mathrm{r})`$ $`=`$ $`{\displaystyle \frac{1}{k_{\mathrm{medium}}}}{\displaystyle \underset{n=1}{\overset{N_{\mathrm{max}}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}p_{nm}𝐍_{nm}^{(1)}(k\mathrm{r})+q_{nm}𝐌_{nm}^{(1)}(k\mathrm{r}),`$ (20)
$`𝐇_{\mathrm{int}}(\mathrm{r})`$ $`=`$ $`{\displaystyle \frac{1}{k_{\mathrm{particle}}}}{\displaystyle \underset{n=1}{\overset{N_{\mathrm{max}}}{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}c_{nm}\mathrm{𝐑𝐠𝐍}_{nm}(k\mathrm{r})+d_{nm}\mathrm{𝐑𝐠𝐌}_{nm}(k\mathrm{r}).`$ (21)
where $`k_{\mathrm{medium}}`$ and $`k_{\mathrm{particle}}`$ are the wavenumbers of the field in the surrounding medium and inside the particle, respectively.

There are $`4N_{\mathrm{max}}(N_{\mathrm{max}}+2)`$ unknown variables—the expansion coefficients $`c_{nm}`$, $`d_{nm}`$, $`p_{nm}`$, and $`q_{nm}`$. Since the fields are vector fields, each point gives multiple equations—four independent equations per point. We generate a grid of $`2N_{\mathrm{max}}(N_{\mathrm{max}}+2)`$ points with equal angular spacings in each of the $`\theta `$ and $`\varphi `$ directions, giving $`8N_{\mathrm{max}}(N_{\mathrm{max}}+2)`$ independent equations. Equal angle spaced points are used for simplicity, although points uniformly distributed about a sphere would be better .

The values of the VSWFs at these points on the particle surface are calculated, and used in the column-by-column calculation of the T-matrix.

The computation time (which is independent of the particle shape, depending only on the containing radius $`r_0`$) is shown in table 1. The calculations were carried out in MATLAB  on a 1.5 GHz PC. The times taken are reasonable in comparison to computation times for EBCM for particles with no symmetry .

Results of a sample calculation are shown in figure 2, where the diagonal scattering matrix elements are shown, calculated using the PMM T-matrix. The scattering matrix elements $`S_{11}`$ and $`S_{22}`$ are shown for scattering in two different planes; the effect of non-axisymmetry is evident.

The accuracy and validity of the PMM-calculated T-matrix will be essentially the same as the accuracy and validity of the point-matching algorithm used in the calculation. Thus, a detailed analysis of our simple proof-of-principle implementation serves little purpose. It is obviously useful to use the best, sufficiently fast, point-matching code available. In view of the mathematical similarity between the T-matrix method and the point-matching method, it should be a simple task to adapt any PMM code to the task of generating T-matrix columns.

## 4 Conclusions

The point-matching method is suitable for the calculation of the T-matrix for particles with no symmetry, provided that the particles are not too large. The method has the advantage of being extremely simple to implement within a general T-matrix package, since most of the required routines will be shared with the existing T-matrix code. This results from the mathematical formalisms of the T-matrix method and the point-matching method being essentially the same. Any point-matching algorithm can be used, with multiple expansion origins, automatic convergence checks, and so on. Since the PMM uses the same field expansions as the EBCM, the same numerical difficulties are to be expected for scatterers with large aspect ratios; in such cases, multiple expansion origin algorithms will be necessary. The accuracy of the PMM T-matrix will be the same as the PMM which is used to calculate it. Naturally, the usual conditions of applicability of the PMM, such as the validity of the Rayleigh hypothesis, need to be considered.

The PMM explicitly depends on the Rayleigh hypothesis—the assumption that the fields can be represented by the expansions (13)—(15) over all space rather than just outside and inside spherical surfaces circumscribing and inscribing the surface of the scatterer. The validity of this assumption for arbitrary scatterers is unknown. However, the use of an overdetermined system of equation may well extend the method somewhat beyond the strict range of applicability of the Rayleigh hypothesis by providing a least squares approximation of the fields between the circumscribing and inscribing surfaces where the VSWF expansions might be non-convergent. One advantage of relying on the Rayleigh hypothesis is that the fields are given everywhere, including the fields internal to the scatterer (a T<sup>(int)</sup>-matrix can be used to relate the internal and incident fields). This applies generally to methods that make use of the Rayleigh hypothesis, such as the generalised separation of variables method . In contrast to this, the EBCM, which avoids the Rayleigh hypothesis, gives the tangential surface fields on the surface of the scatterer, rather than the internal fields.

The point-matching method lacks the generality of DDA and FDTD. In this respect, the recent discrete dipole moment method T-matrix calculations by Mackowski  are particularly promising.

Lastly, we note again that FDTD may prove to be a useful method for T-matrix calculation since it can be used to calculate T-matrices simultaneously for a range of particle sizes.


text:
# Intensity Interferometry with Anyons

## Abstract

A pairwise correlation function in relative momentum space is discussed as a tool to characterize the properties of an incoherent source of non-interacting Abelian anyons. This is analogous to the Hanbury–Brown Twiss effect for particles with fractional statistics in two dimensions. In particular, using a flux tube model for anyons, the effects of source shape and quantum statistics on a two-particle correlation function are examined. Such a tool may prove useful in the context of quantum computing and other experimental applications where studying anyon sources are of interest.

Anyons are particles that exhibit fractional quantum statistics in two dimensions. Existing somewhere between bosons and fermions, anyons have been applied theoretically to a variety of problems such as the fractional quantum Hall effect (FQHE) Halperin (1984); D. Arovas and Wilczek (1984), high temperature superconductivity Laughlin (1988), supersymmetry Plyushchay (1997); Rausch de Traubenberg and Slupinski (1997), and fault-tolerant quantum computing Kitaev (2003). Experimentally, the quasi-particle excitations seen in the FQHE have been shown to possess anyon properties Clark et al. (1988); Simmons et al. (1989).

Fractional spin in two dimensions directly addresses the topological interpretation of quantum statistics. For identical particles in three or more dimensions, the appropriate homotopy group is the permutation group Wu (1984) and there are only two homotopy classes. These correspond directly to fermions and bosons.

However, in two dimensions, the corresponding homotopy group is the braid group Wu (1984) and an infinite number of homotopy classes are possible. As particles are exchanged in relative position, each successive winding of the particles around each other cannot be smoothly deformed into a finite number of configurations as is the case in higher ($`>`$2) dimensions. Each winding corresponds to a separate homotopy class which in turn can be categorized as a different classification of quantum statistics. The term anyon was coined by Wilczek Wilczek (1982a, b) to describe objects having such topological properties, a physical example of which is a charged particle bound to a magnetic flux tube in two dimensions.

This relationship between topology and statistics was glimpsed by Aharonov and Bohm Aharonov and Bohm (1959) and later realized, formalized, and expanded by others Leinaas and Myrheim (1977); G.A. Goldin and Sharp (1980, 1981); Wilczek (1982a, b); Forte (1992); Khare (1997). Because the homotopy group is the braid group, this makes many multi-particle ($`>`$2) anyon problems nearly intractable. An essential analytical difficulty arises because, unlike the case with the permutation group, the multi-particle wave functions for anyons cannot in general be written in a simple way in terms of the single-particle wave functions. Multi-particle anyon wave functions are, in effect, permanently entangled. However, it is exactly this topological property of anyons that makes them robust against decoherence and thus desirable candidates for qubits in fault-tolerant quantum computing Kitaev (2003).

Because quantum interference effects can be sensitive to quantum statistics, it is natural to inquire about the role of anyons in this context. As mentioned, the Aharonov-Bohm effect was an early probe into this fascinating problem. More recently, first order interference effects in a Mach-Zender-style interferometer for non-Abelian anyons in the presence of Aharonov-Bohm flux sources have been derived from applications of the braid group Overbosch and Bais (2001).

In this work, intensity interferometry in momentum space, a second-order interference effect, is explored for incoherently emitted, non-interacting Abelian anyons. This is conceptually related to the recent q-Bose gas interferometry approaches for pions in heavy ion collisions Anchishkin et al. (2000a, b). However, q-bosons and q-fermions are not restricted to two dimensions. In that case, the generalization of quantum statistics is achieved somewhat differently than for anyons Haldane (1991).

Intensity interferometry is a useful tool for characterizing the source geometry of incoherently emitted particles. In addition, the interference effect, expressed in terms of a second-order correlation function, might be regarded as an entanglement measure. Such a tool may prove useful in the context of quantum computing and other experimental applications where studying anyon sources are of interest.

Intensity interferometry was originally developed by Hanbury–Brown and Twiss (HBT) as an alternative to Michelson interferometry to measure the angular sizes of stars in radio astronomy Hanbury Brown and Twiss (1956). By correlating intensities rather than amplitudes, the measurement is insensitive to high frequency fluctuations that would normally make Michelson interferometry prohibitive. When applied to classical waves, the HBT effect is essentially a beat phenomenon.

The field of modern quantum optics was spawned when intensity interferometry was quantum mechanically applied to photons rather than classical waves. Now, two- and multi-photon effects are routinely studied in photonics. Two-fermion HBT in 2D condensed matter systems has also been reported Henny et al. (1999). The HBT effect was independently applied to pions and other particles in elementary particle physics in momentum space and is often called the Goldhaber-Goldhaber-Lee-Pais (GGLP) effect in that context Goldhaber et al. (1960); Boal et al. (1990); Baym (1998).

Quantum mechanically, the two-particle correlation function is a measure of the degree of independence between a joint measurement of two particles ($`1`$ and $`2`$ below). Formally, the correlation function can be written:

$$C_2=\frac{\mathrm{Tr}[\rho a^{}a^{}aa]}{\mathrm{Tr}[\rho a^{}a]\mathrm{Tr}[\rho a^{}a]}\frac{P(1,2)}{P(1)P(2)}$$
(1)
where $`\rho `$ is the density matrix and $`a^{}`$ and $`a`$ are the creation and annihilation operators for the quanta associated with the appropriate fields of interest. When normalized to the single particle distributions as shown, $`C_2`$ is proportional to the relative probability for a joint two-particle measurement as compared to two single-particle measurements. If the measurements are independent, then $`C_2`$=1. If the measurements are correlated, $`C_2`$ deviates from unity. In principle, $`C_2`$ can be expressed as a joint measurement as a function of any degree of freedom of interest. In this treatment, we will focus on momentum space correlations.

The second-order correlation function in momentum space is a powerful tool to probe several important properties of a system. $`C_2`$ is sensitive to the quantum statistics of the particles as expressed in the commutation relations for the field operators in Eq. (1). Also, as determined by the form of the density matrix, $`C_2`$ contains information about the space-time distribution of the particle source in phase space as well as the pairwise interaction and the quantum field configuration.

Non-interacting identical particles can exhibit strong correlations, no correlations, or even anti-correlations depending on the specific field configuration. Possible field configurations for identical bosons include thermal states, coherent states, Fock states, and squeezed states Scully and Zubairy (1997). For example, an incoherent thermal field configuration gives $`C_2=2`$ for spinless bosons of equal momentum. In contrast, a coherent source of bosons, such as laser light well above the lasing threshold, gives a constant $`C_2=1`$ for all relative momenta.

Geometric information about the source is contained in the shape of the correlation function as a function of the pair’s relative momentum. The correlation function typically approaches unity for large momentum differences. For incoherent sources, the relationship between source width and correlation width is essentially limited by the Heisenberg uncertainty principle such that $`\delta q1/R`$ where $`\delta q`$ is the width of $`C_2`$ in relative momentum space and $`1/R`$ is some characteristic width of the incoherent source generating the pairs.

Keeping in mind the above discussion, this paper examines the behavior of $`C_2`$ in momentum space for a non-expanding, non-relativistic, incoherent source of anyons. An approximation to Eq. (1) is used to highlight some features of intensity interferometry in momentum space. The quantum statistics are “tuned”, and, given some simple source functions, the corresponding correlation function is extracted.

In momentum space for an incoherent source, Eq. (1) reduces to the Koonin-Pratt equation Koonin (1977); Pratt et al. (1990); Brown and Danielewicz (1997)

$$C_2(𝐪)1=𝑑𝐫K(𝐪,𝐫)S(𝐫)$$
(2)
where the integration kernel is given by

$$K(𝐪,𝐫)=|\mathrm{\Phi }_𝐪(𝐫)|^21.$$
(3)
The function $`\mathrm{\Phi }_𝐪(𝐫)`$ is the two-particle wave function in the center of mass frame of the pair, where $`𝐪=(𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐})/2`$ is the relative momentum and $`𝐫=𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}`$ is the relative separation in that frame. The center of mass motion of the pair is not considered in this treatment. The source function, $`S(𝐫)`$, is the normalized probability distribution of emitting a particle pair with relative separation $`𝐫`$. The integral is over the entire relative separation space. In this context, incoherent means that the particles are emitted from the source randomly and independently. Moreover, any potential time dependence of the source function in Eq. (2) has been integrated out as in Ref. Brown and Danielewicz (1997). This integration eliminates the ability to recover emission-time ordering information between the pairs and places the focus on spatial information only.

Using a flux tube model for anyons, a dynamical approach is used to obtain two-particle non-relativistic wave functions for anyons in two dimensions Aharonov and Bohm (1959); Wilczek (1982a, b); Khare (1997). In this model, the quantum statistics are enforced through an interaction, the strength of which enforces the effective quantum statistics. In this picture, each anyon is analogous to a charged particle bound to a tightly bundled magnetic field in two-dimensions. The analogy is not strictly one-to-one in this treatment because the scalar Coulomb potential between the charged pairs is not considered; indeed, the system need not be electromagnetic. Nevertheless, the analogy paints a helpful physical picture.

The tightly bundled magnetic field associated with each anyon, $`\stackrel{}{B}=\times \stackrel{}{A}`$, points perpendicular to the plane of motion. We can select a gauge such that the vector potential, $`\stackrel{}{A}`$, has the form $`A_r=0`$ and $`A_\varphi =\xi /2\pi r`$, where $`\xi `$ is the magnetic flux through the plane. The particles are only permitted in regions of space where the magnetic field of the other particle is zero, so the particles exert no force on each other. Quantum mechanically, the interaction term takes the form of a minimal coupling such that the momentum operator takes the form $`\widehat{p}(ie\stackrel{}{A})`$. While under the conditions described above, this “interaction” exerts no force, it does adjust the phase of the wave functions, permitting interference effects. The relative wave function of the pair can then be used in Eq. (2), in combination with a normalized source distribution, to obtain $`C_2`$.

In the center-of-mass frame, using the Hamiltonian

$$\widehat{H}=\frac{\widehat{p}_r^2}{m}+\frac{(\widehat{p}_\varphi \alpha )^2}{mr^2},$$
(4)
the wave function is obtained for the relative motion of two free anyons in two-dimensional polar coordinates. In Eq. (4), $`m`$ is the mass of one particle ($`m_1=m_2=m`$) and $`\alpha `$ represents the anyon parameter with $`0\alpha 1`$, where $`\alpha =0(1)`$ corresponds to bosons (fermions). In this standard treatment, all of the details of the vector potential, magnetic flux, and gauge choice are included in the tunable parameter $`\alpha `$. In Eq. (4), the reduced mass of the identical pair, $`\mu =m_1m_2/(m_1+m_2)=m/2`$, has already been substituted. In the two-dimensional relative coordinate system of the identical pairs, the x-axis has been chosen along the direction of $`\stackrel{}{q}`$.

Applying Eq. (4) to the time-independent Schrödinger equation, $`\widehat{H}\psi =E\psi `$, gives

$$[(\frac{^2}{r^2}+\frac{1}{r}\frac{}{r})\frac{1}{r^2}(i\frac{}{\varphi }+\alpha )^2+q^2]\psi =0$$
(5)
where $`q^2=mE`$ is the square of the relative momentum of the pair (i.e. the energy is that of a free relative particle) and $`\varphi `$ is the angle between $`𝐪`$ and $`𝐫`$. Eq. (5) admits solutions of the form

$$\psi (r,\varphi )=NJ_{|l\alpha |}(qr)e^{il\varphi }$$
(6)
where $`N`$ is a normalization constant. We want solutions that correspond to the appropriate symmetric (antisymmetric) boson (fermion) free particle wave function in relative coordinates in the limit of $`\alpha =0(1)`$. Specifically,

$$\mathrm{\Phi }_{\alpha =0,1}(r,\varphi )=\frac{1}{\sqrt{2}}(e^{i𝐪𝐫}\pm e^{i𝐪𝐫})$$
(7)
where the $`\alpha =0`$ case corresponds to the symmetric solution (+) while the $`\alpha =1`$ case is the antisymmetric solution (-).

Keeping in mind the identity

$$e^{i𝐪𝐫}=e^{qr\mathrm{cos}\varphi }=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}i^nJ_n(qr)e^{in\varphi }$$
(8)
one can construct partial wave superpositions of Eq. (6) that give Eq. (7) for the limiting cases of $`\alpha =0,1`$. Namely,

$$\mathrm{\Phi }_{q,\alpha }(r,\varphi )=\sqrt{2}\underset{l=\mathrm{}(\mathrm{even})}{\overset{+\mathrm{}}{}}(i)^{|l\alpha |}J_{|l\alpha |}(qr)e^{il\varphi }$$
(9)
where $`l`$ is the orbital angular momentum quantum number. The sum is over partial waves of Bessel functions of the first kind of order $`|l\alpha |`$ where $`l`$ is necessarily an integer, but $`\alpha `$ is not. The subscript $`q`$ has been added to the wave function at this stage to remind the reader that we will be integrating over $`r`$ and $`\varphi `$ using Eqs. (2) and (3) to obtain a function of $`q`$.

Treating the initial pair wave function in Eq. (5) as that of spinless bosons (the “boson basis”) simplifies the problem of obtaining the appropriate exchange symmetry for the final two-particle wave function shown in Eq. (9). In constructing Eq. (9), the partial wave expansion using Eq. (6) is symmetrized under the exchange $`\varphi (\varphi \pi )`$. Tuning the anyon parameter will result in a wave function that possesses the correct exchange symmetry for identical indistinguishable pairs. For the fermion and boson cases, Eq. (9), with the help of Eq. (8), reduces to Eq. (7). But for $`0<\alpha <1`$, Eq. (9) is the non-trivial two-anyon free particle wave function.

Equation (2) is used to obtain an expression of $`C_2`$ for anyons. While intensity interferometry can be used to image three dimensional sources, for simplicity we will only consider angle-averaged sources that are a function of $`r`$ rather than $`𝐫`$. The integration kernel in Eq. (3) is angle-averaged in two dimensions and then denoted $`K_{0_\alpha }(q,r)`$ such that

$$K_{0_\alpha }(q,r)=[\frac{1}{2\pi }_0^{2\pi }𝑑\varphi |\mathrm{\Phi }_{q,\alpha }(r,\varphi )|^2]1.$$
(10)
Inserting the expression for the relative wave function from Eq. (9) into Eq. (10) and performing the angle-averaging in two dimensions gives

$$K_{0_\alpha }(q,r)=[2\underset{l=\mathrm{}(\mathrm{even})}{\overset{l=+\mathrm{}}{}}|J_{|l\alpha |}(qr)|^2]1.$$
(11)
In the limiting case when $`\alpha =(0,1)`$, Eq. (11) simplifies to, with the help of Bessel function properties,

$$K_{0_{\alpha =0,1}}(q,r)=(1)^\alpha J_0(2qr)$$
(12)
which is, by construction, the same angle-averaged kernels one gets by using Eq. (7) directly. For the more general case when $`0<\alpha <1`$, the sum in Eq. (11) is evaluated numerically to an appropriately converging order.

Two sample source functions are used to illustrate the effects of source shape and size on $`C_2`$:

$$S_g(r)=\frac{e^{r^2/4r_0^2}}{4\pi r_0^2}$$
(13)
$$S_s(r)=\frac{\mathrm{\Theta }(r_0r)}{\pi r_0^2}.$$
(14)
They represent incoherent pair emission distributions normalized in two dimensions,

$$2\pi _0^{\mathrm{}}r𝑑rS(r)=1,$$
(15)
and where the effective source width is given by $`r_0`$. The forms of the functions were chosen as examples to represent typical localized incoherent sources.

Integrating Eq. (2) over the azimuthal angle in two dimensions using Eq. (11), and choosing either Eqs. (13) or (14) as the source function, gives an angle-averaged expression for $`C_2`$ parameterized by $`\alpha `$:

$$C_{2_\alpha }(q)1=2\pi _0^{\mathrm{}}r𝑑rK_{0_\alpha }(q,r)S(r).$$
(16)
Figure 1 shows $`C_{2_\alpha }`$ for an incoherent Gaussian source, Eq. (13). Various values of $`\alpha `$ are displayed. The figures are shown plotted versus the quantity $`qr_0`$, which accounts for the shifting momentum scale upon changing the source size. For example, as the width of the source decreases, the correlation function widens and vice versa.

For the boson and fermion cases, the standard HBT results for incoherent sources are recovered. The top curve in Fig. 1 represents the spin zero case and approaches a value of $`C_{2_{\alpha =0}}=2`$ as $`q0`$. Here, it is more likely to measure a pair of non-interacting identical bosons in a state of zero relative momentum than to independently measure each particle with their respective momentum. That is, in the joint measurement, the bosons are not independent. This effect is a function of $`q`$ and is the result of the Bose-Einstein statistics and the symmetric nature of the wave function. Similarly, the lowest curve represents non-interacting indistinguishable fermions and displays an anti-correlation in the joint measurement at low relative momentum, where $`C_{2_{\alpha =1}}=0`$ as $`q0`$. Similar to the boson case, this is also only due to the quantum statistics and indicates again that a measurement of one particle affects the other quantum mechanically. As the relative momentum increases, the value of $`C_{2_\alpha }`$ approaches unity, indicating that the joint measurement becomes less correlated.

Starting from the top boson curve and moving downward are the correlation curves for various values of the anyon parameter. Like fermions, the anyons display an exclusion principle and always anti-correlate in the limit as $`q0`$. However, other than the exclusion at $`q=0`$, there is a natural trend of decreasing correlation as one interpolates between bosons and fermions for values of $`q>0`$. The correlation function also approaches unity for large values of $`q`$.

An example of the sensitivity of $`C_{2_\alpha }`$ to the source shape is shown in Fig. 2. The box profile with a sharp edge introduces a ringing structure into the correlation function.

Experimentally, the measured correlation function provides information about the size and shape of the emission source. Measured correlation functions can be inverted to determine actual source distributions using methods like those discussed in Ref. Brown and Danielewicz (1997). Also, information about anyon pairwise interactions can be extracted by comparing measured correlation functions against predictions made using hypothesized interactions in the Hamiltonian. Moreover, the degree of chaoticity of the source, testing the independent emission assumption, can be extracted by looking at the behavior of the correlation at low relative momentum. Finally, the correlation function itself provides information about the degree of entanglement between pairs versus their relative momentum, a potentially useful tool in the field of quantum computing.

In summary, intensity interferometry in momentum space provides a tool that can be used to experimentally study the properties of anyons and their sources in two dimensions. Further theoretical studies using this framework can be pursued exploring various properties of multi-anyon systems including anyon interactions, non-Abelian anyon correlations, and higher-order multi-anyon correlations.

###### Acknowledgements.

The author would like to thank D.A. Brown, J.L. Klay, S.R. Klein, J. Randrup, R. Vogt, and Z. Yusof for valuable input and stimulating discussions.


text:
# Dijet cross sections in ep collisions: who is afraid of symmetric cuts?

## Abstract

Three widely used scenaria for defining cuts on transverse energies of jets in ep collisions are discussed. All of them are shown to suffer from the same sort of unphysical behaviour when the cut regions are subject to additional constraints. This feature is inherent in the very way dijet cross sections are defined and cannot be avoided. In particular, the symmetric cut scenario is shown to be equally suitable for the comparison with NLO QCD calculations as the asymmetric or sum-like ones.

In analyzing the data on inclusive jet production in electron-proton collisions some lower cut $`E_c`$ (at HERA typically $`E_c5`$ GeV) on the transverse energy of jets is necessary for experimental as well as theoretical reasons. In parts of the multijet phase space the cut on the transverse energies of two jets with highest $`E_T`$ may, however, get into conflict with the infrared safety of the corresponding cross sections.

In the next-to-leading QCD calculations of single and dijet cross sections there are at most three final state partons, which combine to two or three jets. Ordering the jets according to their transverse energies in $`\gamma ^{}`$p center-of-mass system ($`E_T^{(1)}E_T^{(2)}E_T^{(3)}`$), the corresponding phase space, marked by the upper and lower solid straight lines in Figs. 2-4, is given by the inequality $`E_T^{(1)}/2E_T^{(2)}E_T^{(1)}`$. This constraint, given by pure kinematics, must then be combined with additional conditions that might be imposed on jets for experimental and theoretical reasons. Phrasing the constraint on allowed part of phase space in terms of $`E_T^{(1)},E_T^{(2)}`$ does not, however, imply that only genuine dijet variables must be investigated or only the first two jets taken into account. For instance, the inclusive single jet cross section gets in principle contribution even from the third jet.

The aim of this letter is to argue that there is no way of choosing the cuts on $`E_T^{(1)},E_T^{(2)}`$ that would be free from the sort of unphysical behaviour of cross sections noted in . We shall illustrate our claims on the calculation of cross sections in ep collisions at HERA in the region relevant for the analysis in

$$2Q^280\mathrm{GeV},0.1y0.85,E_c=5\mathrm{GeV},$$
(1)
but they are clearly of more general validity. All calculations reported below were obtained by means DISENT NLO Monte-Carlo program  using CTEQ6M set of parton distribution functions (PDF) of the proton. For dijet events there are three classes of distributions that can be measured:

* separate distributions of transverse energies and pseudorapidities

$$\frac{\mathrm{d}\sigma }{\mathrm{d}E_T^{(1)}},\frac{\mathrm{d}\sigma }{\mathrm{d}E_T^{(2)}},\frac{\mathrm{d}\sigma }{\mathrm{d}\eta ^{(1)}},\frac{\mathrm{d}\sigma }{\mathrm{d}\eta ^{(2)}},$$
(2)
* their sums, corresponding to distributions of what is called in  “trigger jets”,

$$\frac{\mathrm{d}\sigma }{\mathrm{d}E_T}\frac{\mathrm{d}\sigma }{\mathrm{d}E_T^{(1)}}+\frac{\mathrm{d}\sigma }{\mathrm{d}E_T^{(2)}},\frac{\mathrm{d}\sigma }{\mathrm{d}\eta }\frac{\mathrm{d}\sigma }{\mathrm{d}\eta ^{(1)}}+\frac{\mathrm{d}\sigma }{\mathrm{d}\eta ^{(2)}},$$
(3)
* or the distributions of some of the combinations characterizing the whole event, like their mean values $`\overline{E}_T(E_T^{(1)}+E_T^{(2)})/2`$ and $`\overline{\eta }(\eta ^{(1)}+\eta ^{(2)})/2`$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}\overline{E}_T},\frac{\mathrm{d}\sigma }{\mathrm{d}\overline{\eta }},$$
(4)
but one is free to use general combinations as well.

Symmetric cut scenario

In this scenario the same cut is imposed on both jets

$$E_T^{(1)}E_c,E_T^{(2)}E_c.$$
(5)
This is quite appropriate for the comparison of data with lowest QCD calculations because the latter involve binary partonic hard processes which lead (in diparton rest frame) to two partons with the same transverse energies. Its use in NLO QCD calculations of dijet cross sections has been questioned in  observing that it leads to unphysical behaviour of the cross section $`\sigma _{sym}(\mathrm{\Delta })`$ corresponding to the integral over the region defined in (5) supplemented with an additional constraint $`E_T^{(1)}E_c+\mathrm{\Delta }`$, represented graphically by the red dotted curve in the upper right plot of Fig. 1.

The resulting dependence, displayed as red solid curve in the left part of Fig. 1, implies that for $`\mathrm{\Delta }`$ approaching zero, $`\sigma _{sym}(\mathrm{\Delta })`$ decreases despite the fact that the corresponding phase space increases.

The source of this unphysical behaviour can be traced back to the fact that at the order $`\alpha \alpha _s^2`$ perturbative QCD does not yield a well-defined result for the double differential cross section $`\mathrm{d}\sigma /\mathrm{d}E_T^{(1)}\mathrm{d}E_T^{(2)}`$ when $`E_T^{(1)}=E_T^{(2)}`$. This comes from the fact that at this order there are mass divergencies of the virtual corrections that must be cancelled by those of the real emissions. For a quantity to be calculable in perturbative QCD one must thus allow for integration over sufficient part of the three parton phase space, where $`E_T^{(2)}<E_T^{(1)}`$.

To understand better the origin of the behaviour of $`\sigma _{sym}(\mathrm{\Delta })`$ we plot in Fig. 2 the distributions of $`E_T^{(1)},E_T^{(2)}`$ and $`\overline{E}_T`$ for $`E_c=5`$ GeV. As expected, the distributions of both $`E_T^{(1)}`$ and $`\overline{E}_T`$ turn negative for $`E_T^{(1)}`$ or $`\overline{E}_T`$ close to $`E_c=5`$ GeV because of the small volume (represented in left parts of Figs. 2-4 by the length of perpendicular solid or oblique dotted lines) of phase space available for the real parton emission. In fact both these distributions diverge to $`\mathrm{}`$ as $`E_T^{(1)}`$ or $`\overline{E}_T`$ approach $`E_c`$ from above, thereby causing $`\sigma _{sym}(\mathrm{\Delta })`$ and $`\sigma _{sum}(\mathrm{\Delta })`$ in Fig. 1 to bend down for small $`\mathrm{\Delta }`$. On the other hand $`\mathrm{d}\sigma /\mathrm{d}E_T^{(2)}`$ grows monotonously with decreasing $`\mathrm{\Delta }`$, reflecting the fact that even for $`E_T^{(2)}=E_c`$ there is plenty of phase space (represented in the left parts of Figs. 2-4 by the dash horizontal lines) available for the real gluon emission because $`E_T^{(1)}`$ can be anywhere above it. Note that the integrals over all three distributions in Fig. 2 are finite and equal as they correspond to three different ways of performing the integral over the same region defined in (5). There is thus no reason to reject the symmetric cut scenario for the comparison of data with NLO QCD, but for $`E_T`$ close to $`E_c`$, one should compare the distributions of $`E_T^{(2)}`$, rather than those of $`E_T^{(1)}`$.

Asymmetric cut scenario

To avoid the above unphysical behaviour of $`\sigma _{sym}`$ the asymmetric cut scenario

$$E_T^{(1)}E_c+\delta ,E_T^{(2)}E_c,$$
(6)
where $`\delta 0`$, was suggested in  111Instead of (6) slightly different modification of (5) was used even earlier in ..

In Fig. 3 we plot the corresponding results for the same three distributions as in Fig. 2. As expected, the distribution of $`E_T^{(1)}`$ now behaves physically down to $`E_c+\delta `$, as do those of $`E_T^{(2)}`$ and $`\overline{E}_T`$. However, whereas the first distribution ends there, the other two continue below $`E_c`$. Both of them are, however, singular at $`E_T^{(2)}=\overline{E}_T=E_c`$, diverging at this point from left to $`+\mathrm{}`$ due to the divergence of real emission contributions as $`E_T^{(2)}E_T^{(1)}`$. However, as in the case of the symmetric cuts, the integrals over the distributions in Fig. 3 are finite and the same. The behaviour of the distributions of $`E_T^{(2)}`$ and $`\overline{E}_T`$ below $`E_c`$ is of course as unphysical as the negative values of the distributions of $`E_T^{(1)}`$ and $`\overline{E}_T`$ closely above $`E_c`$ in the symmetric cut scenario of Fig. 2.

Moreover, it is simple to invent an additional constraint that cuts off part of the phase space region defined in (6) in a similar way as the condition $`E_T^{(1)}E_c+\mathrm{\Delta }`$ in the case of the symmetric cut. For instance the restriction

$$E_T^{(2)}2E_T^{(1)}E_c\delta \mathrm{\Delta },$$
(7)
which cuts off part of the region (6) as shown in the middle right picture in Fig. 1 results in the integrated cross section $`\sigma _{asym}(\mathrm{\Delta })`$ exhibiting the same kind of unphysical behaviour for $`\mathrm{\Delta }0`$ as $`\sigma _{sym}(\mathrm{\Delta })`$.

Sum cut scenario

The cuts on $`E_T^{(1)},E_T^{(2)}`$ can also be formulated in terms of their mean value

$$\overline{E}_T\frac{E_T^{(1)}+E_T^{(2)}}{2}E_c+\frac{\delta }{2},E_T^{(2)}E_c.$$
(8)
As shown in Fig. 4 in this case it is the distributions in $`E_T^{(1)}`$ and $`E_T^{(2)}`$ that are singular at $`E_c`$: the former turns negative and diverges when $`E_T^{(1)}`$ approaches $`E_c`$ from above, whereas the latter diverges to plus infinity for $`E_T^{(2)}E_c`$ from below. As in the previous two cases the integrals over all three distributions are, however, finite and the same.

Although even for the lowest $`\overline{E}_T`$ the region defined in (8) involves integration over finite neighbourhood of the point $`E_T^{(1)}=E_T^{(2)}`$, this scenario was claimed in  not to be infrared safe. The argument was based on the fact that similarly as in the case of the symmetric cut scenario (5) the integrated cross section $`\sigma _{sum}(\mathrm{\Delta })`$, corresponding to the additional constraint $`E_T^{(1)}E_c+\delta /2+\mathrm{\Delta }`$ behaves, as illustrated by the green dashed curve in Fig. 1, unphysically for $`\mathrm{\Delta }0`$. However, as emphasized above, the same objection can actually be raised against any scenario. This is clear from the geometrical meaning of the additional cuts involving the parameter $`\mathrm{\Delta }`$: as shown in the right part of Fig. 1, in all three scenaria it cuts off the peaks of the wedges in the plane $`E_T^{(1)},E_T^{(2)}`$ where the negative infinities of the virtual corrections dominate.

The scenario (8) is often supplemented with a constraint on the relative difference of jet transverse energies, for instance, $`(E_T^{(1)}E_T^{(2)})/(E_T^{(1)}+E_T^{(2)})1/4`$. This cuts off part of the three jet phase space but as indicated by the black dash-dotted curve in Fig. 1, has only small effect of the behaviour of $`\sigma _{sum}(\mathrm{\Delta })`$.

The results presented in Figs. 1-4 demonstrate that there is no real advantage of the asymmetric or sum-like cut scenaria over the symmetric one. All of them suffer from the same kind of unphysical behaviour of two of the three distributions in $`E_T^{(1)},E_T^{(2)}`$ or $`\overline{E}_T`$ and all three can be supplemented by additional constraints that lead to the same unphysical dependence on integrated cross sections $`\sigma _{sym}(\mathrm{\Delta })`$, $`\sigma _{asym}(\mathrm{\Delta })`$, or $`\sigma _{sum}(\mathrm{\Delta })`$ displayed in Fig. 1.

Distributions in bins of finite width

In practice experimental distributions in $`E_T^{(1)},E_T^{(2)}`$ and $`\overline{E}_T`$ have finite bin widths. For each bin of these variables we can introduce an additional cut $`\mathrm{\Delta }`$, which chops off part of this bin, and investigate the dependence of the bin content on this parameter in the limit $`\mathrm{\Delta }0`$, much in the same way as we did above for the integrated cross sections. Note that for $`E_T^{(1)},E_T^{(2)}E_c+\delta `$ the choice of the cut scenario is irrelevant for these considerations as it does not influence the distributions there. Let us take as an example the bin

$`7E_T^{(1)}10\mathrm{GeV},`$ $`E_T^{(2)}E_c`$ (9)
$`7E_T^{(2)}10\mathrm{GeV},`$ $`E_T^{(1)}E_T^{(2)}`$ (10)
represented in left two pictures of Figs. 5,6 by the vertical and horizontal colored bands. We first introduce the parameter $`\mathrm{\Delta }`$ similarly as for the symmetric cut (6) by imposing additional constraint:

$$7+\mathrm{\Delta }E_T^{(1)}10\mathrm{GeV}$$
(11)
which means graphically cutting off the red areas in Fig. 5. The resulting $`\mathrm{\Delta }`$-dependence of the contributions $`\mathrm{d}\sigma /dE_T^{(1)}`$ and $`\mathrm{d}\sigma /dE_T^{(2)}`$ to the corresponding bin content $`\sigma _{bin}(\mathrm{\Delta })`$, displayed in the right part of Fig. 5, shows that the distribution of $`E_T^{(2)}`$ behaves unphysically for small $`\mathrm{\Delta }`$, whereas that of $`E_T^{(1)}`$ increases monotonously as $`\mathrm{\Delta }0`$.

However, we can cut off part of the bin (9-10) in an alternative way represented by the red areas in Fig. 6

$$7E_T^{(2)}10\mathrm{\Delta }\mathrm{GeV},$$
(12)
which is certainly as legitimate as the cut (11). The corresponding results, shown in the right part of Fig. 6, lead, however, to opposite conclusion! Now it is the $`\mathrm{d}\sigma /dE_T^{(2)}(\mathrm{\Delta })`$ distribution that is a monotonously decreasing function of $`\mathrm{\Delta }`$, whereas $`\mathrm{d}\sigma /dE_T^{(1)}(\mathrm{\Delta })`$ exhibits much the same unphysical dependence on $`\mathrm{\Delta }`$ as did $`\mathrm{d}\sigma /dE_T^{(2)}(\mathrm{\Delta })`$ for the cut defined in (11).

In summary, for dijet sample there is thus no way how to introduce the cuts on the jet transverse energies that would avoid entirely the problem noted in . Whatever variables, quantity and binning we choose, there is always a simple way how to cut the chosen kinematical region in a way that leads to the same sort of unphysical behaviour of cross sections shown in Fig. 1 and the asymmetric cut scenario is therefore not superior to the symmetric one. On the other hand the above exercise also suggest a simple cure: one just needs to choose the appropriate variable or take sufficiently wide bins in jet transverse energies.

We are grateful to G. Grindhammer and V. Shekelyan for raising the questions that initiated this investigation. This work was supported in part by the project LN00A006 of the Ministry of Education of the Czech Republic and by the Institutional research project AV0Z1-010-920.