text:
# Planets in Stellar Clusters Extensive Search. I. Discovery of 47 Low-amplitude Variables in the Metal-rich Cluster NGC 6791 with Millimagnitude Image Subtraction Photometry.

## 1 Introduction

Since antiquity mankind has wondered whether planetary systems other than our own exist. The first documented effort aimed at extrasolar planet detection was undertaken by Huygens in the XVIIth century. Starting in the 1930s, subsequent searches have been attempted, but failed to produce positive results due to insufficient measurement precision. Only recently it has become possible to obtain radial velocity measurements accurate enough to indirectly detect planets via Doppler shifts of stellar spectra of stars other than the Sun (Mayor & Queloz 1995).

To date, over 70 planets have been discovered, mainly around solar type stars. All of them were detected by radial velocity surveys (eg. Marcy et al. 2001, Noyes et al. 1997). For one of these systems, HD 209458, the transit of the planet across the host star’s disk has been observed (Charbonneau et al. 2000, Henry et al. 2000), thus demonstrating the feasibility of detecting planets this way. Several groups are currently monitoring the brightness of thousands of stars to search for planets via transits (eg. Brown & Charbonneau 1999, Quirrenbach et al. 1998). Recently, Udalski et al. (2002) discovered 46 stars with transiting low-luminosity companions. These objects may be planets, brown dwarfs or M dwarfs. A confirmation of their nature will be provided by mass determinations based on photometry combined with radial velocities from followup spectroscopic observations.

The analysis of the properties of stars with planets suggests that they are on the average significantly more metal rich than those without (Santos et al. 2001). Some studies indicate that the source of the metallicity is most likely “primordial” (Santos et al. 2001, Pinsonneault et al. 2001), while others suggest that the observed high metallicity is intrinsic only in some cases, with the more likely cause being the accretion of planetesimals onto the star (Murray & Chaboyer 2001) or the infall of giant gas planets (Lin 1997).

Although the observed lack of planets in the low metallicity (\[Fe/H\]$`=0.7`$) globular cluster 47 Tuc (Gilliland et al. 2000) is compatible with the “primordial” metallicity scenario, the case is far from being resolved. In such dense environments as the cores of globular clusters encounters with other stars may lead to the breakup of planetary systems (Davies & Sigurdsson 2001). The study of open clusters offers the possibility of disentangling the effects of metallicity and crowding: their stellar densities are not high enough for crowding-induced disruption to be effective.

We have undertaken a long-term project, Planets in Stellar Clusters Extensive Search (PISCES), to search for transiting planets in open clusters. As our first target we have chosen NGC 6791 $`[(\alpha ,\delta )_{2000}=(19^h20.8^m,+37^{}51^{})]`$, a very old, extremely metal rich cluster ($`\tau `$=8 Gyr, \[Fe/H\]=+0.4; Chaboyer et al. 1999). At a distance modulus of (m-M)V = 13.42 (Chaboyer et al. 1999) it contains about 10000 stars (Kaluzny & Udalski 1992, hereafter KU92).

In this paper we present the results of a test observing run at the FLWO 1.2 m telescope. Our aim was to demonstrate the feasibility of obtaining the required accuracy using image subtraction photometry on data collected with a 1-m class telescope needed to reliably detect transits of inner-orbit gas-giant planets with an acceptably low false alarm rate.

On clear nights we have reached the desired level of photometric precision. Unfortunately, there were few such nights during our 26 night run, which is typical for July at FLWO. Even though this dataset is not optimal for planetary transit detection, it allowed us to discover 47 new low amplitude variables, compared to 22 previously known (Kaluzny & Rucinski 1993; Rucinski, Kaluzny & Hilditch 1996, hereafter: RKH).

The paper is arranged as follows: §2 describes the observations, §3 summarizes the reduction procedure, §4 outlines the procedure of variable star selection, §5 contains the variable star catalog. Concluding remarks are found in §6.

## 2 Observations

The data analyzed in this paper were obtained at the FLWO 1.2 m telescope using the 4Shooter CCD mosaic with four thinned, back side illuminated AR coated Loral $`2048^2`$ CCDs (Szentgyorgyi et al. 2002). The camera, with a pixel scale of $`0\stackrel{}{\mathrm{.}}333`$ pixel-1, gives a field of view of $`11\stackrel{}{\mathrm{.}}4\times 11\stackrel{}{\mathrm{.}}4`$ for each chip. The cluster was centered on chip 3 (Fig. 1). The data were collected during 18 nights, from July 6th to August 1st, 2001. None of the nights were photometric and most were at least partially cloudy, as this was the monsoon season. A total of $`204\times 900`$ s $`R`$ and $`36\times 450`$ s $`V`$-band exposures were obtained. The median seeing was $`1\stackrel{}{\mathrm{.}}9`$ in R and $`2\stackrel{}{\mathrm{.}}1`$ in V.

## 3 Data Reduction

### 3.1 Image Subtraction Photometry

The preliminary processing of the CCD frames was performed with the standard routines in the IRAF ccdproc package.111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the NSF.

The photometry for all stars in the field was extracted using the ISIS image subtraction package (Alard & Lupton 1998, Alard 2000) from the $`R`$ and $`V`$-band data for all four CCD chips.

ISIS is based on the fast optimal image subtraction (OIS) algorithm. In order to successfully subtract two images, it is necessary to exactly match their seeing. In OIS this is accomplished by finding a convolution kernel $`Ker`$ and difference in background levels $`bg`$ which will minimize the squared differences between both sides of the equation:

$`Im(x,y)=Ker(x,y;u,v)Ref(u,v)+bg(x,y)`$
where $`Ref`$ is the reference image and $`Im`$ is the processed image. (Alard & Lupton 1998, Wozniak 2000).

The ISIS reduction procedure consists of the following steps: (1) transformation of all frames to a common $`(x,y)`$ coordinate grid; (2) construction of a reference image from several best exposures; (3) subtraction of each frame from the reference image; (4) selection of stars to be photometered; (5) extraction of profile photometry from the subtracted images.

All computations were performed with the frames internally subdivided into four sections ($`sub\mathrm{\_}x=sub\mathrm{\_}y=2`$). Differential brightness variations of the background were fit with a first degree polynomial ($`deg\mathrm{\_}bg=1`$). A convolution kernel varying quadratically with position was used ($`deg\mathrm{\_}spatial=2`$).

An image of particularly good quality was selected as the template frame for the stellar positions. The remaining images were re-mapped to the template frame coordinate system using a second degree polynomial transform (degree=2). During this step an initial rejection of cosmic rays was also performed. A setting of 1.0 for the cosmic ray threshold (cosmic\_thresh) was used.

A reference image was then constructed from 20 best images in $`R`$ and 15 in $`V`$. The constituent images were processed to match the template PSF and background level and then stacked by taking a median in each pixel to obtain a reference image virtually free of cosmic rays.

Image subtraction was then applied to all the frames. For each frame the reference image was convolved with a kernel to match its PSF and then the frame was subtracted from it. As the flux of the non-variable stars on both images should be almost identical, such objects will disappear from the subtracted image. The only remaining signal will come from variable stars.

Next comes the selection of stars for photometry. Normally at this stage the standard ISIS variable detection procedure is used. It computes a median of absolute deviations on all subtracted images and performs a simple rejection of cosmic rays and defects. This approach works very well for variables with an amplitude of at least 0.1 mag, but is much less efficient for variables of smaller amplitude. As we are interested in these latter variables we extracted photometry for all the stars on the template list, to search them for variability using more traditional methods.

A final step, not included in the ISIS image subtraction package itself, was to convert the light curves from ADU to magnitudes. For this purpose the template frames were reduced with the DAOphot/Allstar package (see Section 3.2 for details). The template instrumental magnitude $`m_{tpl}`$ of each star was converted into counts $`c_{tpl}`$, using the Allstar zeropoint of 25 mag. The light curve was then converted point by point to magnitudes $`m_i`$ by computing the total flux $`c_i`$ for a given epoch as the sum of the counts on the template $`c_{tpl}`$ and the counts on the subtracted template image $`\mathrm{\Delta }c_{tpl}=c_{ref}c_{tpl}`$ decreased by the counts corresponding to the subtracted image $`\mathrm{\Delta }c_i=c_{ref}c_i`$:

$$c_i=c_{tpl}+\mathrm{\Delta }c_{tpl}\mathrm{\Delta }c_i$$
(1)
The flux was converted to instrumental magnitudes using the same zero point as above. To convert the photometric error $`\sigma _i^c`$ expressed in counts to $`\sigma _i^m`$ in magnitudes, we used the following relation:

$$\sigma _i^m=2.5\mathrm{log}\left(\frac{c_i}{c_i+\sigma _i^c}\right)$$
(2)
### 3.2 Profile Photometry

Profile photometry was extracted using the DAOphot/Allstar package (Stetson 1987) from images remapped onto the template $`(x,y)`$ coordinate grid. This was done mainly to compare its precision to that of image subtraction photometry.

A PSF varying linearly with position on the frame was used. Stars were identified using the FIND subroutine and aperture photometry was performed with the PHOT subroutine. The same star list was used on all frames for the construction of the PSF. Of those the stars with profile errors greater than twice the average were rejected and the PSF was recomputed. This procedure was repeated until no such stars were left on the list. The PSF was then further refined on frames with all but the PSF stars subtracted from them. This procedure was applied twice. The resultant PSF was then used by Allstar in profile photometry.

The next step was to obtain a template list of stars, to be also used as the input star list for the photometry routine in ISIS. The template and the stacked reference image were reduced following the same procedure as for single frames. The stars were then subtracted from the reference image and an additional 100-200 stars undetected in the initial FIND run were hand-selected and added to the star list. Allstar was run again twice on the expanded list, which resulted in the rejection of some more objects. The star list was further cleaned by rejecting objects with $`\chi >3`$ and with errors too large for their magnitude. This list was used as input to Allstar on the template frame to create the final list of stars.

The final template star list was then used as input to Allstar in the fixed-position mode on each of the frames. The output profile photometry was transformed to the common instrumental system of the template image and then combined into a database. The database was created for the $`R`$-band images in chip 3 only.

### 3.3 Implications for transit detection

In Figure 2 we plot the ISIS (left panel) and DAOphot/Allstar (right panel) time series precision as a function of magnitude for the best night. This night should be representative of the average night we expect other than during the July-August monsoon season. The time series precision was computed as the rms of the light curve, with the rejection of $`3\sigma `$ outliers. Only stars with light curves containing at least 22 of 23 data points are plotted. The continuous curve indicates the photometric precision limit due to Poisson noise of the star and sky brightness. The dashed curves correspond to the rms value which would yield a 6.5 $`\sigma `$ detection for three transits in the $`R`$-band, for planets with 1, 1.5, 2 and 3 $`R_J`$ and orbital period of 3.5 days, the same as for HD 209458. The curves are defined by the equation

$$rms=\sqrt{N}\frac{\mathrm{\Delta }R}{\sigma }$$
(3)
where $`N`$ is the number of observations during three transits, $`\mathrm{\Delta }R`$ is the $`R`$-band amplitude of the transit and $`\sigma `$ is set at 6.5. The steep fall of detection capability at the bright end is caused by the increase in radius on the subgiant branch and the 8 hour limit on the transit length set by the length of the night. The increase at the faint end is caused by the decrease of the stellar radius along the main sequence.

Not surprisingly, the overall precision is better and the scatter smaller for image subtraction photometry. This is especially well seen for the brightest stars. In the plot for image subtraction photometry 6562 (95%), 5304 (77%), 4110 (59%) and 2053 (29%) stars fall within the detection limit for transits of 3, 2, 1.5 and 1 $`R_J`$ planets, respectively.

There is a group of stars with $`rms>0.1`$ mag and $`R<20`$, with several discrepant points in the light curve caused by bad columns, which were not removed by sigma clipping. Following the approach of Udalski et al. (2002) we will set an upper limit on the $`rms`$ to remove such stars from the sample searched for planetary transits.

Figure 3 shows the transit length (upper panel) and the probability of occurence of a transit given random inclinations (lower panel) as a function of $`R`$-band magnitude, plotted for periods of 1, 3.5, 5, 10 and 15 days. The transit length is proportional to the period $`P`$ as $`P^{\frac{1}{3}}`$ and the probability of observing a transit given random inclinations, is proportional to $`P^{\frac{2}{3}}`$. Both of these quantities depend on the mass, $`M`$, and radius, $`R`$, of the star as $`M^{\frac{1}{3}}R`$ (Eqs (1) in Gilliland et al. 2000). The duration of the transits is of the order of 1 to 4 hours and it decreases with the magnitude of the star. For periods over $`3.5`$ days the probability of a transit is of the order of 10% for stars near the turnoff and it falls to 2-5% for lower main sequence stars.

Figure 4 shows the detection efficiency (upper panel) and cumulative detection efficiency (lower panel) for transiting planets with periods greater than $`1^d`$ as a function of period, assuming $`8^h`$ of observations per night during a $`30`$ day continuous run. The choice of the minimum period is motivated by the shortest period transiting object ($`P=0.8082^d`$) reported by Udalski et al. (2002), although it should be noted that the shortest period for a planet detected by radial velocity searches is $`3`$ days (HD 83443, Mayor et al. 2000). The detection efficiency for at least two transits drops below 30% at $`P9^d`$ (upper panel, dot-dashed line) and for at least three transits at $`P5^d`$ (solid line). Such temporal coverage would enable us to observe three or more transits for 50% of all transiting planets with periods between $`1`$ and $`7.8`$ days (lower panel, solid line). A 50% efficiency for the detection of at least two transits is reached at a period of $`15.7`$ days (dot-dashed line).

The distribution of best night residuals for 6871 stars, normalized to the rms, is shown in Fig. 5. Only stars with at least 22 out of 23 data points were included in the histogram. The distribution of the residuals is well-fit by a Gaussian distribution with $`\sigma =0.96`$.

### 3.4 Calibration

None of the nights in our observing run were photometric, thus we were forced to rely on other, indirect sources of calibration.

The $`R`$-band photometry in chip 3 was calibrated from the photometry for 12 of the 14 red giant branch (RGB) stars published by von Braun et al. (1998). One of the stars could not be identified on our template and another one, much redder than the others (R4, $`\mathrm{V R}=1.355`$) was rejected because it was a clear outlier in both $`R`$ and $`V`$. Due to the small number of stars and a very limited color baseline, only an offset was determined. The residuals, with an rms scatter of 0.014 mag, are shown in the upper panel of Fig. 6. The $`R`$-band photometry for chips 1, 2 and 4 was left uncalibrated, as they contain no stars with calibrated $`R`$-band photometry.

A similar comparison was made for the $`V`$-band magnitudes of the 12 RGB stars and there seems to be a trend with color (lower panel of Fig. 6) and magnitude, as those two quantities are correlated for RGB stars. A trend in magnitude is also present for stars with $`V<17`$ in comparison with data of Kaluzny & Rucinski (1995; hereafter KR95), but is not seen in case of the KU92 photometry. Since the KR95 data cover parts of all four chips and the offset is in good agreement with the von Braun et al. (1998) photometry (to $`0.005`$ mag), we have used offsets from the KR95 photometry to bring the instrumental magnitudes to the standard system.

In Figures 7 and 8 we present the $`V/\mathrm{V R}`$ CMDs for chip 3 and chips 1, 2, 4, respectively. The CMD for chip 3, centered on the cluster, shows a well defined main sequence (MS) down to $`V20`$, a subgiant branch (SGB), a red giant branch and a red clump at $`V14.5`$, $`\mathrm{V R}0.7`$. The blue stars, noted by KU92, are also present at $`\mathrm{V R}<0`$. The CMDs for the remaining chips consist in large proportion of disk stars; only an indication of the upper MS is discernible.

### 3.5 Astrometry

Equatorial coordinates were determined for the $`R`$-band template star lists, expanded with the variables with no $`R`$-band photometry. The transformation from rectangular to equatorial coordinates was derived using 734, 729, 1058 and 657 transformation stars from the USNO A-2 catalog (Monet et al. 1996) in chips 1 through 4, respectively. The average difference between the catalog and the computed coordinates for the transformation stars was $`0\stackrel{}{\mathrm{.}}14`$ in right ascension and $`0\stackrel{}{\mathrm{.}}12`$ in declination.

## 4 Variability Search

### 4.1 Selection of Variables

To select candidate variable stars we have employed the variability index J defined by Stetson (1996). Below we present only a summary of the method; the reader is referred to the original paper for details.

The variability index J is computed as follows:

$$J=\frac{_{k=1}^nw_k\mathrm{sgn}(P_k)\sqrt{|P_k|}}{_{k=1}^nw_k}$$

where the user has defined n pairs of observations, i and j, to be considered, each with a weight $`w_k`$. In our case, observations separated by less than 0.03 days were treated as a pair. A weight $`w_k=1`$ was assigned to pairs of observations ($`i(k)j(k)`$) and $`w_k=0.25`$ to single observations ($`i(k)=j(k)`$). $`P_k`$ is the product of the normalized residuals of the two observations, i and j, constituting the k-th pair:

$$P_k=\{\begin{array}{cc}\delta _{i(k)}\delta _{j(k)},\hfill & \mathrm{if}i(k)j(k)\hfill \\ \delta _{i(k)}^21,\hfill & \mathrm{if}i(k)=j(k)\hfill \end{array}$$

and $`\delta `$ is the magnitude residual of a given observation from the average, scaled by the standard error:

$$\delta =\sqrt{\frac{n}{n1}}\frac{v\overline{v}}{\sigma _v}$$

The final variability index was multiplied by a factor $`w/w_{max}`$, where $`w_{max}`$ is the total weight of a star if it were measured on all images.

We are aware that this approach is not optimal for detecting planetary transits and a better method based on the matched-filter algorithm is under development. Recently, Udalski et al. (2002) used this technique to identify 46 transiting planet candidates. They report that their implementation was very sensitive even to single transit events and produced virtually no spurious detections.

We believe the choice of the algorithm was not critical for the dataset analized here. Due to the poor weather conditions during most of the observing run and the resulting decreased photometric accuracy on some nights and uneven temporal coverage, we did not expect to find planetary transits.

### 4.2 Period Determination

To search for periodicities we have used the method introduced by Schwarzenberg-Czerny (1996), employing periodic orthogonal polynomials to fit the observations and the analysis of variance statistic to evaluate the quality of the fit.

If observations $`X`$ consist of the sum of signal $`F`$ and noise $`E`$: $`X_k=F_k+E_k`$, then the analysis of variance statistic $`\mathrm{\Theta }(\omega )`$ is defined as:

$$\mathrm{\Theta }(\omega )\frac{\widehat{Var}\{F\}}{\widehat{Var}\{E\}}=\frac{(K2N1)_{n=0}^{2N}|c_n^2|}{2N[(X,X)_{n=0}^{2N}|c_n^2|]}$$

where 2N is the order of the complex polynomial (corresponding to a Fourier series of N harmonics), K is the number of observations and $`c_n`$ are the coefficients of the orthonormal polynomial $`\mathrm{\Psi }_N(z)=_{n=0}^Nc_n\mathrm{\Phi }_n(z)`$, where $`\mathrm{\Phi }_N(z)=_{n=0}^Na_n^{(N)}z^n`$ is the base. For the details of the method the reader is referred to the Schwarzenberg-Czerny (1996) paper.

## 5 Variable Star catalog

We have found 62 variables in all four chips, 47 of them new ones (24 in chip 3): 30 eclipsing binaries, 21 other periodic variables and 11 miscellaneous ones. Their $`R`$-band light curves are shown in Figures 9, 10 and 11 and parameters listed in Tables 1, 2 and 3, respectively.222The $`R`$ band photometry and finding charts for all variables are available from the authors via the anonymous ftp on cfa-ftp.harvard.edu, in the /pub/bmochejs/PISCES directory. The variables are also plotted on CMDs for their respective chips \[Figs 7 (chip 3) and 8 (chips 1, 2, 4)\]. Eclipsing binaries are marked with circles, other periodic variables with squares and miscellaneous variables with triangles.

We have recovered 18 of the 22 previously known variables. The missing four are V7 and V8, which are outside of our feld of view, and V13 and V19, which are saturated in the $`R`$-band data. The light curves of the long period detached eclipsing binaries V10, V18 and V21 are also not shown, as we have not observed any eclipses for them.

We have reobserved the W UMa type eclipsing variables V1-V6. We confirm the longer of the two possible periods found for V5 by RKH. This variable shows unequal maxima, as noted previously by RKH. The shape of the light curve is seen to have changed: in our data the variable attains higher maximum brightness at phase 0.25 (between the primary and the secondary eclipse) and in the RKH data at phase 0.75. The light curve of V4 shows a similar change: the maxima are almost equal now, while in the RKH data they were clearly uneven. The variability of the maxima is probably due to changing starspots on the surface of one or both of the components, like in RS CVn type binaries. We have discovered six new W UMa type systems in chips 1, 2 and 4, V22-V28. One of them, V23, also shows unequal maxima.

V6 exhibits extremely small amplitudes of 0.1 mag in $`R`$ and $`V`$, considering its flat-bottomed minima. Such light curve shape would indicate an inclination angle close to $`90^{}`$ and an extremely small mass ratio for the binary components ($`q<0.05`$, Rucinski 1993). A comparison with the 0.15 mag amplitude light curve in Kaluzny et al. (1993) strongly suggests that this variable is blended in our data. This is confirmed by the elongated appearance of the star on our reference frame.

Another short period eclipsing variable, V29, does not seem to be a contact system judging by the narrowness of the primary eclipse. On the CMD it is located below the main sequence, together with V2 and V11. It is most likely a field star behind the cluster.

B4, one of the blue stars identified by KU92, shows brightness variations with an amplitude of 0.05 mag. We have phased the light curve with a period of 0.8 days, so that it displays two minima and two maxima. B4 might be a non-eclipsing binary with brightness modulations due to ellipsoidal variations or a reflection effect. An argument in favor of this interpretation comes from a study of the radial velocities of a sample of 70 sdBs which indicates that 45% of them are post-common envelope binaries with periods of the order of a few hours to a few days (Saffer, Green & Bowers 2000). The period for B4 is consistent with this period range.

The binary nature of B4 is interesting from the point of view of the origin of hot subdwarfs in metal rich clusters. Binary evolution has been suggested as one of the possible formation mechanisms of such objects (Moehler 2001).

We have also found some variables which could either be EB type eclipsing binaries or ellipsoidal variables, as most have amplitudes smaller than 0.1 mag. Among them are V14 and V16 with newly derived periods. The new variables which belong to this type are V30, V31, V33, V34, V36 and possibly V37 and V38. During one of the nights V37 seems to have undergone a UV Ceti-like outburst with an amplitude of 0.1 mag (Fig 9). This variable is also very red, with V R=1.1 and most likely not a cluster member.

We have reobserved the RS CVn type eclipsing binary V9. The shape of the S-wave does not seem to have changed significantly from the one observed by RKH in 1996.

We have derived periods for the EA type eclipsing binaries V11, V12. For the detached eclipsing system V20 we have observed only one eclipse. Three new such sytems, V35, V60 and V61, were observed as well. A period for V35 was derived, as it appears that one of the eclipses was observed twice. Other possible EA type binaries are V32 and V39.

In addition to eclipsing binaries and ellipsoidal variables we have also discovered periodic variables which seem to be of different nature. The newly phased variable V17, as well as the new variables V40-V44, V46, V47, V51, V53 are most likely BY Dra type variables – rotating spotted late type dwarfs. On the CMDs they are located in the vicinity of the lower MS ($`V>18.5`$) or redward of it. Some are probably not cluster members. V40-V42 and possibly V46 exhibit small humps between the maxima. Similar behavior has been observed for other BY Dra type variables, eg. CC Eri (Amado et al. 2000). We have also considered the possibility that V40-V42 are field RR Lyrae stars. Despite the similarity in the light curve shape and period, this interpretation seems unlikely because of their small amplitudes.

The remaining periodic variables are brighter than $`V=18`$. The nature of their variability is unclear. These variables are spread throughout the CMDs. V48, V52 and V58 seem to be located at the cluster turnoff, and V45, V54, V56, V49 and V50 above it. V55, significantly redder than the RGB stars, is probably not a cluster member. The variables V57 and V59 seem to be located redward of the turnoff. V59 exibits quite a high amplitude of 0.3 mag.

Among the miscellaneous variables there are two especially interesting ones: the cataclysmic variables B7 (V15) and B8 (Kaluzny et al. 1997). B7 was observed to undergo a 0.6 mag outburst while in its high state. B8 showed two dwarf-novae type outbursts, the first one with an amplitude of 2 mag and duration of about 5 days and a second, slightly fainter one.

Five of the six other miscellaneous variables, V63-V67, are located above the cluster turnoff. V62 appears to be located on the MS. Most of the bright variables seem to show long period or quasi-periodic variability. V66 and V67, located $`0.1`$ mag redward of the RGB, display variations on a longer timescale than the other variables.

## 6 Conclusions

In this paper we have demonstrated the feasibility of obtaining photometry accurate enough to detect planets through transits in open clusters with 1 m class telescopes. The analysis of the data collected for this purpose resulted in the discovery of 47 new low amplitude variables, compared to 22 previously known (Kaluzny & Rucinski 1993, RKH).

The first stage of our project will be to obtain for this and several other clusters continuous observations with a 1-m class telescope for about $``$30 nights under very good observing conditions. The photometry will be obtained in two bands, $`R`$ and $`V`$, to allow us to differentiate between planetary transits and blended eclipsing binaries. Planetary transits should be gray, as the planet does not contribute any measurable light to the system, while the superposition of an eclipsing binary with another star will show a change in color during the eclipse.

In the next stage better light curves will be obtained with a 2 m class telescope for selected planet transit candidates. Using radial velocity measurements derived from spectroscopic observations with a 6-10 m class telescope it will be possible to distinguish planetary and brown dwarf transits from grazing eclipses by main sequence companions. A precision of 1 km s-1 will enable us to identify and reject stars with companions above 0.0075 M ($`8`$ MJ).

We would like to thank Janusz Kaluzny for helpful discussions and careful reading of the manuscript and Alceste Bonanos for her help in obtaining some of the data. We also thank the referee, Ronald Gilliland, for a very useful report which significantly improved the paper. Part of this work was done during BJM’s visit to the Harvard-Smithsonian Center for Astrophysics. BJM was supported by the Polish KBN grant 5P03D004.21 and the Foundation for Polish Science stipend for young scientists. DDS acknowledges support from the Alfred P. Sloan Foundation and from NSF grant No. AST-9970812.


format: text

text:
# submitted to The Astrophysical Journal Letters Do the Infrared Emission Features Need Ultraviolet Excitation? The PAH Model in UV-poor Reflection Nebulae

## 1 Introduction

Since their first detection in the planetary nebulae NGC 7027 and BD+30o3639 (Gillett, Forrest, & Merrill 1973), the distinctive set of infrared (IR) emission features at 3.3, 6.2, 7.7, 8.6, and 11.3$`\mathrm{\mu m}`$ have been observed in a variety of objects with a wide range of physical conditions, including planetary nebulae, protoplanetary nebulae, reflection nebulae, HII and ultracompact HII regions, circumstellar envelopes, the diffuse interstellar medium (ISM) of the Milky Way Galaxy, and external galaxies (see Tielens et al. 1999 and Sellgren 2001 for recent reviews). Although the exact nature of their carriers remains unidentified – they remain known as “the Unidentified Infrared (UIR) bands” – it is now widely thought that they originate from some sorts of aromatic hydrocarbons. Various carbonaceous materials have been proposed as the UIR band carriers. In general, they can be categorized into two classes: (1) pure, free-flying aromatic molecules – polycyclic aromatic hydrocarbon molecules (PAHs; Léger & Puget 1984; Allamandola, Tielens, & Barker 1985); (2) carbonaceous grains with a partly aromatic character – hydrogenated amorphous carbon (HAC; Duley & Williams 1981; Jones, Duley, & Williams 1990), quenched carbonaceous composite (QCC; Sakata et al. 1990), coal (Papoular et al. 1993), fullerenes (Webster 1993), and surface-graphitized nanodiamonds (Jones & d’Hendecourt 2000).

In most current models, the UIR band emission involves three sequential steps: (1) excitation by absorption of an energetic starlight photon (usually ultraviolet \[UV\]); (2) rapid ($`10^{12}10^{10}\mathrm{s}`$) redistribution of all or part of the absorbed photon energy over all available vibrational modes; (3) radiative relaxation via IR fluorescence.

Among the existing proposed carriers, the PAH model is gaining increasing acceptance because of (1) the close resemblance of the UIR spectra (frequencies and relative intensities) to the vibrational spectra of PAH molecules (e.g. see Allamandola, Hudgins, & Sandford 1999) and (2) the ability of a PAH molecule to emit efficiently in the UIR wavelength range following single photon heating (Greenberg 1968; Léger & Puget 1984; Allamandola, Tielens, & Barker 1985; Draine & Li 2001). In contrast, larger carbonaceous grains are unlikely to be efficient UIR emitters since the timescale for the absorbed photon energy to diffuse ($`10^9\mathrm{s}`$) is much shorter than the IR vibrational emission timescale ($`0.1\mathrm{s}`$; see Tielens et al. 1999).

Early observations of the UIR emission bands were made in regions with strong UV irradiation, and PAH excitation was expected since all PAH species are strongly absorbing in the vacuum ultraviolet ($`\lambda 3000`$Å). The PAH UV excitation model has recently been challenged by the ISO (Infrared Space Observatory) detection of the UIR bands in $`\mathrm{vdB}\mathrm{\hspace{0.17em}133}`$, a reflection nebula illuminated by a binary system with little UV radiation (see Uchida, Sellgren, & Werner 1998). The UIR spectrum of this UV-poor region closely resembles those observed in sources with much harsher UV environments (Uchida et al. 2000). This appears to be in conflict with the view that PAHs are primarily excited by UV photons, as would be expected based on laboratory studies showing that the absorption by small, neutral PAHs has a sharp cut-off in the UV, with little or no absorption in the visible (see Sellgren 2001 for a review). Uchida et al. (2000) obtained 5–15$`\mathrm{\mu m}`$ ISOCAM spectra of 6 reflection nebulae and found no systematic spectroscopic differences despite values of $`T_{\mathrm{eff}}`$ ranging from 6800$`\mathrm{K}`$ to 19000$`\mathrm{K}`$.

The PAH electronic absorption edge is known to shift to longer wavelength with increasing size and/or upon ionization (see Allamandola et al. 1989, Salama et al. 1996 and references therein). While the largest experimentally studied PAH molecule to date is dicoronylene C48H20 (Allamandola, Hudgins, & Sandford 1999), astronomical PAHs are believed to be larger (e.g., the mean size for the Milky Way PAHs is $`6`$Å, corresponding to $`N_\mathrm{C}100`$ \[Li & Draine 2001a\]). The astronomical PAH model – with the size/ionization dependence of the PAH absorption edge taken into account (see §A2 in Li & Draine 2001a) – is successful in explaining the observed mid-IR spectra of the Milky Way diffuse ISM (Li & Draine 2001a), the quiescent molecular cloud SMC B1#1 in the Small Magellanic Cloud (Li & Draine 2001b), and the UIR band ratios for a wide range of environments ranging from reflection nebulae, HII regions, photodissociation regions (PDRs), molecular clouds in the Milky Galaxy to normal galaxies, starburst galaxies, and a Seyfert 2 galaxy (Draine & Li 2001).

In this Letter we show that the astronomical PAH model is consistent with the observed UIR emission from UV-poor environments. In §2 we verify that the astronomical PAH model can quantitatively reproduce the vdB 133 UIR spectrum. We further demonstrate in §3 that the UIR bands are also expected for reflection nebulae which are even more UV-poor than vdB 133, and we provide model spectra for comparison with future Space Infrared Telescope Facility (SIRTF) observations of reflection nebulae near cool stars. In §4 we show that the predicted $`T_{\mathrm{eff}}`$ dependence of the ratio of the IRAS 12$`\mathrm{\mu m}`$ emission to the total far-IR emission is also consistent with observations. Our results are discussed in §5, and our conclusions are summarized in §6.

## 2 vdB 133

The reflection nebula vdB 133, located in the Cyg OB1 association at a distance $`d1500\mathrm{pc}`$, is illuminated by a binary system consisting of HD 195593A (spectral type F5Iab; effective temperature $`T_{\mathrm{eff}}6800\mathrm{K}`$; absolute magnitude $`M_V=6.6\pm 0.4`$) and HD 195593B (spectral type B7II; $`T_{\mathrm{eff}}12000\mathrm{K}`$; $`M_V=4.0\pm 0.4`$) (see Uchida et al. 2000 and references therein). Following Uchida et al. (1998), we approximate the vdB 133 radiation field by a linear combination of two Kurucz model atmosphere spectra (Kurucz 1979). The starlight intensity incident on the reflection nebula is assumed to have $`U_0135`$, where $`U_0`$ is the ratio of the 912Å–1$`\mathrm{\mu m}`$ energy density relative to the value for the MMP interstellar radiation field ($`4.80\times 10^{13}\mathrm{erg}\mathrm{cm}^3`$).

Let $`f_\lambda `$ be the fraction of the nonionizing stellar luminosity $`L_\lambda `$ radiated at wavelength $`<\lambda `$:

$$f_\lambda \frac{_{912\mathrm{\AA }}^\lambda L_\lambda ^{}𝑑\lambda ^{}}{_{912\mathrm{\AA }}^{\mathrm{}}L_\lambda ^{}𝑑\lambda ^{}}$$
(1)
This binary system,111HD 195593B, with $`f_{0.2\mu \mathrm{m}}32\%`$, contributes $`80\%`$ of the total flux shortward of 0.2$`\mathrm{\mu m}`$; HD 195593A, with $`f_{0.2\mu \mathrm{m}}0.4\%`$, contributes $`95\%`$ of the total flux shortward of 1$`\mathrm{\mu m}`$.  with $`f_{0.2\mu \mathrm{m}}1.7\%`$, has the lowest value of $`f_{0.2\mu \mathrm{m}}`$ for reflection nebulae with detected UIR emission bands.222 For comparison, NGC 2023 is illuminated by HD 37903, with $`T_{\mathrm{eff}}=22000\mathrm{K}`$ and $`f_{0.2\mu \mathrm{m}}68\%`$; the average interstellar starlight background (Mathis, Mezger, & Panagia 1983) has $`f_{0.2\mu \mathrm{m}}=4.7\%`$.  Despite this, the vdB 133 UIR features are very similar to those observed in regions with much harder radiation fields (Uchida et al. 2000).

We model the vdB 133 UIR features as emitted from a mixture of neutral and ionized PAH molecules. We use the PAH absorption cross sections estimated by Li & Draine (2001a) which take into account the size/charge dependence of the PAH long wavelength absorption edge (see Figure 1). We adopt a log-normal size distribution for the PAHs, characterized by three parameters: $`a_0`$, $`\sigma `$, and $`b_\mathrm{C}`$; $`a_0`$ and $`\sigma `$ respectively determine the peak location and the width of the log-normal distribution, and $`b_\mathrm{C}`$ is the total amount of C atoms (relative to total H) locked up in PAHs. As always, the term “PAH size” refers to the radius $`a`$ of a hypothetical spherical grain with the same carbon density as graphite ($`2.24\mathrm{g}\mathrm{cm}^3`$) and containing the same number of carbon atoms $`N_\mathrm{C}`$: $`a1.288N_\mathrm{C}^{1/3}\mathrm{\AA }`$.

The PAH ionization fraction depends on the radiation intensity, the electron density $`n_e`$, and the gas temperature $`T_{\mathrm{gas}}`$. Assuming $`T_{\mathrm{gas}}=100\mathrm{K}`$, we estimate the PAH ionization fractions in vdB 133 using the rates for photoelectric emission and electron capture recently discussed by Weingartner & Draine (2001).

We employ the “thermal-discrete” method (Draine & Li 2001) to calculate the temperature distribution functions for PAHs subject to single-photon heating. Our model assumes a single radiation field: attenuation of the illuminating starlight in the dust layer is neglected. If the dust layer is optically-thick to the exciting radiation, we may approximate this by an optically-thin layer with a column density $`N_\mathrm{H}`$ such that the optical depth to the exciting radiation would be unity. For a standard dust mixture and $`\lambda =4000\mathrm{\AA }`$ (say), this would correspond to $`N_\mathrm{H}1.3\times 10^{21}\mathrm{cm}^2`$.

Figure 2 plots the theoretical spectra calculated from models with $`a_0=3.5`$Å, $`\sigma =0.4`$, and abundances $`(N_\mathrm{H}/10^{21}\mathrm{cm}^2)b_\mathrm{C}=38,33,29\mathrm{ppm}`$ for $`U_0/n_e=250,500,1000\mathrm{cm}^3`$ respectively. We see that the $`U_0/n_e=500\mathrm{cm}^3`$ model ($`n_e=0.3\mathrm{cm}^3`$) closely reproduces the observed UIR spectrum using the standard Milky Way PAH size distribution ($`a_0=3.5`$Å, $`\sigma =0.4`$), a realistic radiation field, and PAH abundances ($`b_\mathrm{C}33\mathrm{ppm}`$ if $`N_\mathrm{H}=10^{21}\mathrm{cm}^2`$) comparable to the Milky Way abundances ($`b_\mathrm{C}45\mathrm{ppm}`$).

## 3 Cool Reflection Nebulae

We have seen in §2 that the astronomical PAH model is successful in explaining the UIR spectrum observed in vdB 133. In this section we will show that reflection nebulae illuminated by cooler stars are also capable of exciting PAHs sufficiently to emit at the 6.2, 7.7, 8.6 and 11.3$`\mathrm{\mu m}`$ UIR bands. For illustrative purpose, we consider reflection nebulae illuminated by stars with $`T_{\mathrm{eff}}=3000,3500,4500,5000\mathrm{K}`$. The latter 3 cases will serve as a direct comparison basis for SIRTF observations of vdB 135 ($`T_{\mathrm{eff}}=3600\mathrm{K}`$), vdB 47 ($`T_{\mathrm{eff}}=4500\mathrm{K}`$), and vdB 101 ($`T_{\mathrm{eff}}=5000\mathrm{K}`$) (Houck 2001).333Attempts by Uchida et al. (2000) to search for UIR bands in vdB 101 and vdB 135 using ISOCAM spectroscopy were unsuccessful due to the faintness of these nebulae.

Lacking prior knowledge of the PAH size distribution and the physical environmental parameters such as $`n_e`$, $`T_{\mathrm{gas}}`$, and the starlight intensity $`U_0`$, we adopt the Milky Way PAH size distribution (i.e., a log-normal distribution with $`a_0=3.5`$Å and $`\sigma =0.4`$ \[Li & Draine 2001\]) and we take $`T_{\mathrm{gas}}=100\mathrm{K}`$. We take $`U_0=100`$ as a representative number since $`30U_0300`$ for typical reflection nebulae (see Sellgren et al. 1990); the calculated UIR emission is essentially proportional to $`U_0`$ for $`U_010^4`$ (see Figure 3a). In all models we take the number of carbon atoms relative to hydrogen in PAHs to be $`b_\mathrm{C}=45\mathrm{ppm}`$, the same as in the Milky Way diffuse ISM (Li & Draine 2001).

For stars with $`T_{\mathrm{eff}}3500\mathrm{K}`$ we take the Kurucz atmospheric model spectra to approximate the radiation fields. Since the Kurucz model is not available for stars cooler than 3500$`\mathrm{K}`$, we take the model spectrum of Allard et al. (1997) for $`T_{\mathrm{eff}}=3000\mathrm{K}`$.

In Figure 3a we plot the mid-IR spectra calculated for PAHs heated by a $`T_{\mathrm{eff}}=3000\mathrm{K}`$ star with $`U_0=1,100,10^3,10^4`$ and $`U_0/n_e=10^4`$. The predicted PAH IR bands are prominent, with a much stronger 11.3$`\mathrm{\mu m}`$ C-H out-of-plane band in comparison with those in regions illuminated by much hotter stars (e.g., NGC 2023 \[$`T_{\mathrm{eff}}=22000\mathrm{K}`$\], Orion Bar \[$`T_{\mathrm{eff}}=37000\mathrm{K}`$\], M17 \[$`T_{\mathrm{eff}}=45000\mathrm{K}`$\]; Verstraete et al. 2001). This is due to the lack of UV photons capable of causing photoelectric ejection of electrons and thus PAHs in such a cool region mainly remain neutral even for $`U_0/n_e=10^4`$. As expected from the single-photon heating scenario, the PAH IR features are almost indistinguishable for $`U_010^4`$ except the $`U_0=10^4`$ model has a rising wing for $`\lambda >10\mathrm{\mu m}`$. This is because for $`U_0=10^4`$ the $`>20\mathrm{\AA }`$ dust grains already attain an equilibrium temperature $`T100\mathrm{K}`$.

In Figure 3b,c,d we plot the model spectra for $`U_0=100`$, $`n_e=0.01,100\mathrm{cm}^3`$, and $`T_{\mathrm{eff}}=3500,4500,5000\mathrm{K}`$, respectively. It is seen that all spectral shapes are almost the same except the absolute emissivity level (at fixed $`U_0`$) increases with increasing of $`T_{\mathrm{eff}}`$. This is because hotter stars radiate a larger fraction of their power at short wavelengths where PAHs are strongly absorbing. It is also seen from the enhanced 11.3$`\mathrm{\mu m}`$ band relative to the 7.7$`\mathrm{\mu m}`$ C-C stretching band that a large fraction of PAHs in those regions are neutral.444The probability for a small PAH molecule being neutral is larger for $`U_0/n_e=10^4\mathrm{cm}^3`$ than that for $`U_0/n_e=1.0\mathrm{cm}^3`$. This is because in regions devoid of UV photons the charging of a neutral PAH is dominated by electron capture and thus PAHs are more likely to be negatively charged when $`U_0/n_e`$ is small.  This is the major difference between the UIR features excited by cool stars with those excited by hot stars.555The vdB 133 nebula is not lack of photoelectron-ejecting photons. The hot companion HD 195593B ($`T_{\mathrm{eff}}=12000\mathrm{K}`$) of the exciting binary system contributes $`80\%`$ of the total 6–13.6$`\mathrm{eV}`$ photons. This is why the vdB 133 UIR spectrum resembles those excited by hot stars.  It is hoped that SIRTF observations of these reflection nebulae will provide both spectra of the UIR bands and a full spectral energy distributions to allow the illuminating intensity $`U_0`$ to be well-determined, thereby allowing comparison with the model spectra in Figure 3.

In Table 1 we give the PAH mid-IR emission convolved with the 5.8$`\mathrm{\mu m}`$ and 8.0$`\mathrm{\mu m}`$ SIRTF Infrared Array Camera (IRAC) bands calculated for the Milky Way PAH mixture excited by cool stars with $`T_{\mathrm{eff}}=3000,3500,4500,5000\mathrm{K}`$. Figure 3 and Table 1 may provide guidance for the coming SIRTF observation of cool reflection nebulae (Houck 2001). We stress, however, that the actual mid-IR spectra may differ from those in Figure 3 and Table 1 since the PAH size distribution, the PAH abundance, and the environmental parameters $`n_e`$, $`T_{\mathrm{gas}}`$ and $`U_0`$ may differ from those adopted here.

## 4 The 12$`\mathrm{\mu m}`$ Emission vs. Radiation Hardness

To study the carrier of the 12$`\mathrm{\mu m}`$ emission and its excitation mechanism, Sellgren et al. (1990) investigated how the ratio $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}/I_{\mathrm{bol}}(\mathrm{FIR})`$ depended on the effective temperature $`T_{\mathrm{eff}}`$ of the exciting stars for 24 reflection nebulae, where $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}`$ is estimated from the IRAS 12$`\mathrm{\mu m}`$ filter, and the total far-IR surface brightness $`I_{\mathrm{bol}}(\mathrm{FIR})`$ was estimated by fitting a blackbody multiplied by a $`1/\lambda `$ emissivity to the IRAS 60 and 100$`\mathrm{\mu m}`$ surface brightnesses. They found that $`\lambda I_\lambda (12\mathrm{\mu m})/I_{\mathrm{bol}}(\mathrm{FIR})`$ is basically independent of $`T_{\mathrm{eff}}`$ for $`5000\mathrm{K}T_{\mathrm{eff}}33000\mathrm{K}`$. If PAHs are the 12$`\mathrm{\mu m}`$ carrier and excited only by UV photons, one would expect $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}/I_{\mathrm{bol}}(\mathrm{FIR})`$ to show a sharp drop for cool stars ($`T_{\mathrm{eff}}<10000\mathrm{K}`$) which are relatively devoid of UV radiation (Sellgren et al. 1990). Therefore, one can conclude that either PAHs are not responsible for the 12$`\mathrm{\mu m}`$ emission or PAHs are excited by photons of a wide range of wavelengths. Since it has been well established that the UIR features are the dominant contributor to the 12$`\mathrm{\mu m}`$ emission, the former hypothesis can be ruled out. Therefore, it must be possible for photons longward of the UV to excite the UIR carrier. Is this consistent with our understanding of PAHs?

To test the astronomical PAH model, we calculate $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}/I_{\mathrm{bol}}(\mathrm{FIR})`$ for reflection nebulae illuminated by stars with $`T_{\mathrm{eff}}`$=3000, 3500, 4500, 6000, 8000, 10000, 15000, 22000 and 30000$`\mathrm{K}`$ and starlight intensity $`U_0=1,100,1000`$, assuming $`U_0/n_e=100\mathrm{cm}^3`$ and using the Milky Way dust mixture. As shown in Figure 4a the model results are in good agreement with observational data except four nebulae: vdB 135 ($`T_{\mathrm{eff}}=3600\mathrm{K}`$), vdB 42 ($`T_{\mathrm{eff}}=4200\mathrm{K}`$), vdB 47 ($`T_{\mathrm{eff}}=4500\mathrm{K}`$), vdB 35 ($`T_{\mathrm{eff}}=4900\mathrm{K}`$). It is worth noting that these nebulae are the four faintest among Sellgren et al.’s sample of 24. It is possible that $`I_{\mathrm{bol}}(\mathrm{FIR})`$ was underestimated for these cool, faint nebulae as a consequence of the method of obtaining $`I_{\mathrm{bol}}(\mathrm{FIR})`$ because the 60 and 100$`\mathrm{\mu m}`$ emission preferentially account the warm grain component while the bulk of the far-IR emission is probably longward of 100$`\mathrm{\mu m}`$ for cool reflection nebulae.

A better way to characterize the fraction of total IR emission emitted at 12$`\mathrm{\mu m}`$ is to replace $`I_{\mathrm{bol}}(\mathrm{FIR})`$ by $`I_{\mathrm{bol}}(\mathrm{tot})_0^{\mathrm{}}\mathrm{I}_\lambda d\lambda `$. In Figure 4b we show the model predicted $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}/I_{\mathrm{bol}}(\mathrm{tot})`$ as a function of $`T_{\mathrm{eff}}`$.

## 5 Discussion

It is not surprising that the astronomical PAH model is able to reproduce the vdB 133 UIR spectrum. It is also not unexpected that the PAH mid-IR emission spectra in regions devoid of UV radiation resemble those in UV-rich regions. This can be readily seen from Figure 7 of Draine & Li (2001) which shows that a 0.5$`\mathrm{\mu m}`$ photon will excite a $`N_\mathrm{C}=100`$ (mean size for the Milky Way mixture) PAH molecule to a vibrational temperature of $`T390\mathrm{K}`$, resulting in strong emission in the 6.2–11.3$`\mathrm{\mu m}`$ UIR bands. For smaller PAHs, softer photons are able to fulfill this excitation task. Therefore, the detection of UIR bands in UV-poor regions does not contradict PAHs as the UIR band carrier. The lack of a precipitous drop for cool stars ($`T_{\mathrm{eff}}<10000\mathrm{K}`$) in $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}/I_{\mathrm{bol}}(\mathrm{FIR})`$ simply indicates that interstellar PAHs have absorption extending into the visible, presumably as the combined result of somewhat larger sizes than those currently studied in laboratory, plus modification of the absorption properties by ionization.

Papoular (2000) proposed that PAH excitation could result from energy released during PAH-catalyzed H2 formation. While this excitation mechanism may be present, it does not appear able to provide sufficient excitation to account for the observed UIR emission (Pagani et al. 1999). In any case, our model calculations show that excitation by starlight can account for the observed UIR intensities.

We note that our approach to estimating the absorption properties of astronomical PAHs was simplified (Li & Draine 2001) in that we did not distinguish negatively-charged anions (PAH-) and multiply-charged cations (PAHn+) and anions (PAHn-) from singly-charged cations (PAH+). For example, laboratory studies of the electronic absorption spectra of coronene (C24H12) in solution have shown that its di-anion (C24H$`{}_{}{}^{2}{}_{12}{}^{}`$) has markedly enhanced visual/near-IR continuum absorption than its mono-anion (C24H$`{}_{}{}^{}{}_{12}{}^{}`$) and neutral counterpart (Hoijtink 1959). However, the fact that our simplified PAH model can successfully reproduce the vdB 133 UIR spectrum and the $`T_{\mathrm{eff}}`$ dependence of $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}/I_{\mathrm{bol}}(\mathrm{FIR})`$ suggests that the model provides a fairly good approximation to the real absorption properties of astronomical PAH mixtures.

## 6 Conclusion

We have modeled the excitation of PAH molecules in UV-poor regions. It is shown that the astronomical PAH model provides a satisfactory fit to the UIR spectrum of vdB 133, a reflection nebulae with the lowest ratio of UV to total radiation among reflection nebulae with UIR bands detected. It is also shown that astronomical PAHs can be pumped by cool stars with even less UV radiation. It is further shown that the PAH model predicts a dependence of $`[\lambda I_\lambda (12\mathrm{\mu m})]_{\mathrm{IRAS}}/I_{\mathrm{bol}}(\mathrm{FIR})`$ for reflection nebulae which is consistent with observations for $`3000\mathrm{K}T_{\mathrm{eff}}30000\mathrm{K}`$. We conclude that PAHs appear able to account for the UIR band emission observed in reflection nebulae.

We thank K.I. Uchida for providing us with the ISO spectrum of vdB 133, F. Allard, C. Dominik and P.H. Hauschildt for their help in obtaining stellar model spectra, and R.H. Lupton for the availability of the SM plotting package. We thank L.J. Allamandola, J.M. Greenberg, K. Sellgren and K.I. Uchida for helpful discussions. This research was supported in part by NASA grant NAG5-10811 and NSF grant AST-9988126.


text:
# Transverse Phase Locking for Vortex Motion in Square and Triangular Pinning Arrays

## I Introduction

The phenomena of phase locking can occur in systems interacting with a periodic potential when an external AC drive is superimposed over a DC drive. A resonance or frequency matching can occur between the AC frequency and the internal frequency generated by motion over the periodic substrate. The best known example is that of Shapiro steps in AC/DC-driven single Josephson-junctions  and arrays of Josephson-junctions , where steps are observed in the current-voltage characteristics. Shapiro like phase locking has also been extensively studied in charge-density-wave systems , spin-density waves  and vortices in superconductors with periodic substrates . In addition, in elastic media moving over random disorder, a washboard signal can be generated by the periodicity of the elastic media itself , allowing Shapiro like steps to be observed when an AC drive is applied, such as in vortex systems with random disorder . In all these systems the AC drive is superimposed in the same direction as the DC drive. Another characteristic of Shapiro like phase locking steps is that the step width oscillates with AC amplitude and frequency between zero and a finite value as an $`n`$th order Bessel function $`J_n(A)`$, where $`n=0`$ corresponds to the depinning or critical current and $`n>0`$ corresponds to higher order steps.

Recently a different kind of phase locking was proposed for a 2D system of vortices moving through a periodic substrate. Here the AC drive is applied in the perpendicular direction, transverse to the longitudinal applied DC drive . Phase locking steps, distinct from Shapiro steps, were observed in the velocity versus applied DC drive. The most pronounced step occurs for DC drive values that allow a sinusoidal orbit to fit into each square pinning plaquette. For larger AC amplitudes a number of additional fractional steps were also observed. In addition, unlike Shapiro steps, the width of the depinning and higher order steps were shown analytically and in simulation to increase monotonically with AC amplitude as $`(A/\omega )^2`$, where $`A`$ is the AC amplitude and $`\omega `$ is the AC frequency. A similar transverse phase locking phenomena can also occur for an ordered vortex lattice moving over random disorder, as was recently demonstrated . In the random disorder case, the steps differ from those observed for square periodic pinning. The step width for random disorder is more like a Shapiro step, going as $`|J_1(A)|`$. In addition, for low AC amplitudes, the step width grows linearly, rather than quadratically as in the case of square pinning.

There are several open questions for the transverse phase locking in the square pinning arrays. It is not known how the additional fractional steps scale with AC amplitude, what the vortex orbits look like along these steps, or what the noise signatures are. It also not known if there is any scaling at the transition to the phase locked state similar to that found in charge density wave systems, where critical scaling is known to occur at depinning and along the transitions into and out of the phase locked regions . The effect of disorder on the transverse phase locking has also not been investigated. The vortex filling fraction or applied field in experiments will likely determine the step width or whether the steps occur at all. For certain filling fractions the overall vortex lattice can be disordered or frustrated , which may destroy the steps. In addition typical transport experiments are often done near $`T_c`$ so the effects of finite temperature on the fractional and main steps should also be investigated to determine the feasibility of observing these states experimentally.

It would also be interesting to examine the effects of a transverse AC drive on systems in which there is an internal transverse oscillation even in the absence of an AC drive. An example of such a system is vortex motion through a triangular pinning array. Interstitial vortices driven through such an array by a longitudinal DC drive oscillate back and forth in the transverse direction as they move in order to pass through each minima of the interstitial potential. In this case the additional internal transverse oscillation may lock with the AC transverse force, leading to an enhancement of the transverse phase locking similar to that observed for the square case. Alternatively, the phase locking may be more Shapiro like, such as the steps found in the transverse phase locking with random pinning .

We focus on the transverse phase locking for vortex motion in square and triangular pinning arrays in thin film superconductors. Vortex pinning and dynamics in thin samples with periodic pinning have attracted considerable attention, since the effects of pinning can be systematically controlled. These pinning arrays can be created using an array of holes  or magnetic dots . Various pinning geometries have also been created such as square , triangular , rectangular , and Kagome . Pronounced commensurability effects appear as peaks in the critical current at integer matching fields and fractional fields where the vortices can form a symmetrical configuration with the pinning array so that the vortex-vortex interaction energy is reduced. These matching configurations have been imaged directly in experiment  and also observed in simulations . For large pinning sites above the first matching field, where there are more vortices than pinning sites, multiple vortices can be captured by each pin up to a saturation point. Beyond this saturation number, additional vortices will be located between the pinning sites in the interstitial regions. For small pinning sites, the saturation number is one, so that for all fields above the first matching field the vortices are located at the interstitial regions. Direct imaging and simulations have shown that these interstitial vortices can form ordered commensurate arrangements . The interstitial vortices are still pinned even though they do not interact directly with the pinning sites, since a periodic potential is created by the repulsive interaction of the vortices trapped at the pinning sites. The interstitial vortices are typically far less strongly pinned than the vortices at the pinning sites, so that under an applied drive the interstitial vortices move while the vortices at the pinning sites remain immobile. Evidence for the motion of the interstitial vortices has been provided by transport measurements  as well as direct imaging in square pinning arrays, which show that the interstitial vortices move in 1D paths between the pinning sites . In simulations, 1D interstitial vortex flow between pinning sites has also been observed . Recent experiments  and simulations  have demonstrated the appearance of steps in the V-I characteristics for moving interstitial vortices when a DC and AC drive are superimposed. In the previous theoretical work for an applied transverse AC drive, vortex filling fractions between $`1.0`$ and $`2.0`$ were examined, where each pinning site captures only one vortex .

In this work we examine the additional fractional steps for the square pinning array with a transverse AC drive. We find fractional steps both above and below the main step. We show that the width of most of the fractional steps scales as the square of the AC amplitude, as previously shown for the depinning current and the main step . Along the steps, the vortex orbits are stabilized in fixed trajectories, while in the non-step regions, the vortex motion is chaotic in appearance, producing a distinct velocity noise spectra. We also find scaling in the velocity versus applied DC drive at depinning and near the main step, with $`\beta =1/2`$. The steps are most prominent for commensurate vortex filling fractions where the interstitial vortices form a symmetrical ground state. For filling fractions where the vortices are disordered, the phase locking phenomena is absent, and for filling fractions just off of a commensurate filling, the steps have a linear increase with drive due to the fact that certain portions of the sample are not undergoing phase locking. We find that even for small finite temperatures, the fractional steps become completely washed out. The main step is visible up to the melting point of the interstitial vortex lattice. With triangular pinning we find a much more pronounced main step and fractional steps, with the depinning threshold nearly the same as that of the square case. The widths of the main and most of the fractional steps scale linearly with the AC amplitude for triangular pinning.

## II Simulation

We numerically integrate the overdamped equations of motion, which for a single vortex $`i`$ is

$$𝐟_i=𝐟_i^{vv}+𝐟_i^{vp}+𝐟_{DC}+𝐟_{AC}=𝐯_i$$
(1)
The total force acting on vortex $`i`$ is $`𝐟_i`$, $`𝐟_i^{vv}`$ is the total force from all the other vortices, $`𝐟_i^{vp}`$ is the force from the substrate, $`𝐟_{DC}`$ is the force from the DC drive and $`𝐟_{AC}`$ is the force from the AC drive. The simulation is in 2D with periodic boundary conditions in the $`x`$ and $`y`$ directions. For the square pinning case we place an array of $`N_p=N\times N`$ pinning sites in the sample, where each pinning site is an attractive parabolic potential of radius $`r_p`$ and a maximum force $`f_p`$. To create a triangular pinning array, we displace every other row of pinning sites by half a pinning lattice constant. The force from the pinning sites is $`f_i^{vp}=_k^{N_p}(f_p/r_p)|𝐫_i𝐫_k^{(p)}|\mathrm{\Theta }(r_p|𝐫_i𝐫_k^{(p)}|)\widehat{𝐫}_{ik}^{(p)}`$,

where $`\mathrm{\Theta }`$ is the step function, $`𝐫_k^{(p)}`$ is the location of pinning site $`k`$, and $`\widehat{𝐫}_{ik}^{(p)}=(𝐫_i𝐫_k^{(p)})/|𝐫_i𝐫_k^{(p)}|`$. In our system each pin captures only one vortex. The total force from the vortex-vortex interaction is $`𝐟_i^{vv}=_{ji}^{N_v}_iU_v(r_{ij})`$ where we use $`U(r)=ln(r)`$, the Pearl vortex-vortex interaction potential, appropriate for thin film superconductors. The experiments with periodic pinning arrays are in this 2D limit. For computational efficiency we use a summation technique  to evaluate the vortex-vortex interaction. Temperature can also be applied to the system by adding a noise term $`f_i^T`$ to the equation of motion. This noise term has the properties $`<f_i^T(t)>=0.0`$ and $`<f_i^T(t)f_j^T(t^{^{}})>=2\eta k_BT\delta _{ij}\delta (tt^{^{}})`$. Here we set $`\eta =k_B=1`$.

We have considered pinning arrays from $`4\times 4`$ to $`10\times 10`$ pins and observe the same results in each case. The initial vortex configurations are obtained by annealing from a high temperature where the vortices are in a liquid state and cooling to $`T=0`$. The DC driving is in the $`x`$-direction, $`𝐟_{DC}=f_{DC}\widehat{𝐱}`$, and the AC force is in the $`y`$-direction, $`𝐟_{AC}=A\mathrm{sin}(\omega t)\widehat{𝐲}`$. The DC force is increased from zero in increments of 0.0009. We average the vortex velocities $`V=𝐯_i\widehat{𝐱}`$ at each increment for 175000 MD steps, and the resulting DC force versus velocity curve is proportional to the DC voltage-current curve.

## III Fractional Steps for Square Pinning Arrays

In Fig. 1 we show the vortex velocity $`V`$ versus applied DC drive $`f_{DC}`$ for two different AC amplitudes, $`A=0.25`$ (upper curve) and $`A=0.35`$ (lower curve), for a fixed frequency of $`\omega =0.002`$. Here there is a large step, centered near $`f_{DC}=0.32`$ for $`A=0.25`$ and $`f_{DC}=0.35`$

for $`A=0.35`$, which corresponds to a phase where the vortex orbit is commensurate with the pinning, forming a sinusoidal like orbit as studied previously. We label this step the “main step” or 1:1 step. In addition to this step, a number of smaller or fractional steps are visible both above and below the main step. For example, for $`A=0.35`$ a prominent step appears near $`f_{DC}=0.27`$, and a small step is visible near $`f_{DC}=0.46`$. The widths of all the steps, including the depinning threshold, are smaller for the lower AC amplitude. The depinning threshold drops from $`f_{DC}=0.196`$ for $`A=0.35`$ to $`f_{DC}=0.190`$ for $`A=0.25`$, while the width $`\mathrm{\Delta }`$ of the main step drops from $`\mathrm{\Delta }=0.09`$ for $`A=0.35`$ to $`\mathrm{\Delta }=0.04`$ for $`A=0.25`$. In Fig. 2(a) we show a magnified region of $`0.222<f_{DC}<0.241`$, near the depinning threshold, for the $`A=0.35`$ curve in Fig. 1. Here the large number of fractional steps present are clearly visible, with the step widths getting progressively smaller for lower $`f_{DC}`$. In Fig. 2(b) we show a magnified region of the $`A=0.35`$ curve above the main step, from $`0.4<f_{DC}<0.45`$, where a series of fractional steps appear. Again the step widths decrease as $`f_{DC}`$ is lowered toward the main step.

### A Commensurate Vortex Orbits and Noise

On steps at small $`f_{DC}`$ values, a vortex spends more than one AC cycle in a single plaquette, while on steps at larger $`f_{DC}`$, the vortex moves through more than one plaquette in a single AC cycle. In Fig. 3 we illustrate the vortex trajectories (lines) and the vortex positions (dots) for fixed $`f_{DC}`$ drives on and near the steps in the $`V`$ versus $`f_{DC}`$ curve of Fig. 1 for $`A=0.35`$. Fig. 3(a) shows the vortex trajectory for the step near $`f_{DC}=0.235`$, which is most clearly visible in Fig. 2(a). Here, a stable vortex orbit occurs in which the vortex spends two AC cycles in one plaquette before moving over to the next plaquette to repeat the sequence. We observe similar behavior at the smaller steps for $`f_{DC}<0.235`$, where in the stable orbits, the vortex spends $`n`$ or $`n+1/2`$ complete AC cycles in one plaquette. In Fig. 3(b), in the stable orbit for the step at $`f_{DC}=0.243`$ in Fig. 1, the vortex spends $`3/2`$ of an AC cycle in each plaquette. In Fig. 3(c) we show the vortex

trajectory for $`f_{DC}=0.28`$, which corresponds to the second largest step in Fig. 1. Here the vortex goes through one AC cycle per plaquette. In the non-step region of Fig. 3(d), at $`f_{DC}=0.305`$, there is no single stable vortex orbit. Instead the vortex trajectory wanders over time. There is an area close to the occupied pinning sites where the vortex never flows due to the vortex-vortex repulsion. We find similar orbits for the other non-step regions. On the main step, illustrated in Fig. 3(e) at $`f_{DC}=0.36`$, the vortex moves through two plaquettes in a single AC cycle. In Fig. 3(f) the vortex orbit for a step above the main step, $`f_{DC}=0.429`$, has a similar pattern to that seen on the main step, Fig. 3(e), except that the vortex moves through three plaquettes in each complete AC cycle.

The presence of stable orbits can also be inferred

through measurement of the velocity noise, which corresponds experimentally to the voltage noise. We acquire a time series of the velocity in the $`x`$-direction for a fixed $`f_{DC}`$ value. The corresponding noise spectrum is given by $`S(f)=|V(t)e^{i2\pi ft}𝑑t|^2`$. In Fig. 4, we examine the noise spectra for $`A=0.35`$ for a step region, $`f_{DC}=0.36`$, and a non-step region, $`f_{DC}=0.305`$, corresponding to the orbits in Fig. 3(e) and Fig. 3(d). In the step region (lower curve), a strong periodic signal is present at the frequency induced by the AC driving, as indicated by the narrow large spike and the higher harmonics. For the non-step region (upper curve) there is is a broadened peak near the AC driving frequency. The noise power for the lower frequencies is over 4 orders of magnitude higher than for the step region. This lower frequency noise power comes from the vortex orbits drifting over a time scale that is longer than the AC driving time scale. The noise spectra does not have a broad band or $`1/f^\alpha `$ characteristic. We find a similar noise spectra for the other non step regions.

### B Scaling Velocity vs Drive and Fractional Steps

The velocity curves in Fig. 1 show considerable curvature near the depinning threshold and just above each step. Fisher has suggested that depinning can be considered a dynamical critical phenomena  where

$$V=(ff_c)^\beta ,$$
(2)
where $`f`$ is the DC driving force, $`f_c`$ is the driving force at which depinning occurs, and $`V`$ is the velocity of the driven media. In addition, in numerical work on CDW systems a similar scaling is found near the mode locked steps  in the form

$$|VV_q|=|ff_q|^\beta ,$$
(3)
where $`V_q`$ is the velocity along the $`q`$th step and $`f_q`$ is the driving force at which the step ends or begins. For a single particle moving in a 1D periodic substrate, $`\beta =1/2`$. In charge density wave systems with random pinning,

simulations of the depinning and step transitions give $`\beta =0.67`$ in 2D and $`0.83`$ for 3D, while experimental results are consistent with the results for the 3D regime . In our system, one might expect to find scaling with $`\beta =1/2`$ as in the 1D periodic substrate picture since the vortices are commensurate with the periodic pinning. On the other hand, the orbits in Fig. 3 indicate that the vortices have a 2D velocity component and do not flow in strictly 1D channels, so it is not clear if or how the velocities will scale.

In Fig. 5(a) we plot on a log scale $`V=ff_c`$ just above depinning for $`A=0.25`$, with the depinning threshold of $`f_c=0.196`$, and in Fig. 5(b) we plot $`|VV_q|=ff_q`$ above the main step for $`A=0.25`$, with the velocity $`V_q=0.24`$ along the main step, and where the end of the step

is at $`f_q=0.395`$. Here in both cases for low $`f`$ the curves fit well to $`\beta =1/2`$, with a deviation to linear behavior for higher $`f`$. Our velocity resolution was too low to obtain adequate scaling for steps other than the main step. We find similar scaling for the depinning and the main step for all filling fractions at which the interstitial vortices form a symmetrical configuration. For non-symmetrical vortex configurations, the scaling breaks down due to the fact that the depinning process for the non-symmetrical configurations is plastic, rather than elastic as in the case of the symmetrical configurations. Scaling may still occur for the incommensurate systems as found in simulations with random pinning where plastic flow occurs, where scaling in the velocity versus drive was observed  with $`\beta >1.5`$. Very large systems beyond the scope of this work would be needed in order to determine if such scaling occurs for the incommensurate cases. In addition, some of the scaling breaks down for high AC amplitudes, as seen for $`A=0.35`$ in Fig. 1 where there are sharp jumps into and out of certain steps. For example, in Fig. 6 we plot the $`dV/df_{DC}`$ curve for $`A=0.35`$, showing peaks going into and out of the phase locked regions. For $`\beta =1/2`$ scaling, the overall shape of $`dV/df_{DC}`$ should scale with a power $`1/2`$, which is consistent with the curve in Fig. 6.

We next show how the widths of different steps scale with AC amplitude. In earlier work it was shown that the increase in the depinning force and the width of the main step increase as the square of AC amplitude. In Fig. 7 we plot the step widths $`\mathrm{\Delta }`$ vs AC amplitude $`A`$ corresponding to the steps at $`f_{DC}=0.235`$ \[Fig. 3(a)\], $`f_{DC}=0.243`$ \[Fig. 3(b)\], and $`f_{DC}=0.28`$ \[Fig. 3(c)\]. We also show the scaling for the main step, $`f_{DC}=0.36`$ \[Fig. 3(e)\], for a comparison. The widths all fit well to $`A^2`$ with some deviations at higher AC amplitudes, where the widths increase faster than $`A^2`$.

### C Effect of Filling Fraction

In Fig. 8 we show a series of velocity vs $`f_{DC}`$ curves for increasing filling fractions from (bottom to top) $`B/B_\varphi =`$ $`17/16`$, $`5/4`$, $`0.34`$, $`3/2`$, $`0.67`$, $`1.97`$, $`2.0`$, and $`2.04`$. The interstitial vortex lattice forms an ordered symmetrical ground state for $`B/B_\varphi =`$ $`17/16`$, $`5/4`$, $`3/2`$, and $`2.0`$. At these fillings the phase locking steps are the most pronounced, and the velocity curves are almost identical for each one. These are also the same filling fractions where peaks in the critical current are observed in experiments and simulations. Since the interstitial vortices form a symmetrical ground state, the effective interstitial vortex interactions cancel, and the system can be thought of as a single particle moving through a periodic substrate. In the absence of an AC drive, the vortices collectively move in 1D paths under a DC drive. For the filling fractions $`B/B_\varphi =0.34`$ and $`0.67`$, the ground state is disordered and the phase locking is absent. In the absence of an AC drive the vortices move plastically in winding paths. Under an applied AC drive, this winding motion destroys any phase locking.

For filling fractions $`B/B_\varphi =1.97`$ and $`2.04`$, just below and just above the commensurate filling at $`B/B_\varphi =2.0`$, the phase locking still occurs; however, the steps have an additional linear velocity increase superimposed on them. In Fig. 9(a) we illustrate vortex trajectories along the main step for $`B/B_\varphi =2.0`$ and in Fig. 9(b) for $`B/B_\varphi =2.04`$. For $`B/B_\varphi =2.0`$ the vortices move in well defined stable sinusoidal trajectories, while for $`B/B_\varphi =2.04`$ only certain rows of vortices have this ordered motion while the vortices in the other rows have unstable ordering. Only the rows with stable orbits are phase locked and remain at the same velocity along the step for increasing $`f_{DC}`$. The vortices moving through the unstable channels are not phase locked and therefore their velocity increases with $`f_{DC}`$, accounting for the linear velocity increase in the $`V`$ \- $`f_{DC}`$ curves along the step. The stable moving vortex rows have a commensurate number of vortices. In the system presented here, for an $`8\times 8`$ pinning lattice there are eight interstitial vortices in each commensurate row. There are more than eight vortices in the rows with unstable vortex trajectories, so the rows are incommensurate. For certain filling fractions such as $`B/B_\varphi =0.67`$, all the rows are incommensurate and there is no phase locking. For $`B/B_\varphi =1.97`$, where the unstable or incommensurate rows have less than eight vortices, we observe vortex trajectories similar to those in Fig. 9(b). We note that similar behavior occurs just above and below the filling fractions of $`5/4`$ and $`3/2`$. These results are consistent with experimental results in which Shapiro steps were most clearly defined at $`B/B_\varphi =2.0`$ .

### D Effect of Temperature

We consider the effects of temperature by adding a noise term $`f_i^T`$ to the equation of motion for the vortices. We normalize our temperature by the melting temperature $`T_m`$, where $`T_m`$ is the temperature at which the vortices began to diffuse in the absence of an external drive. In our system this onset is sharp and well defined, due

to the fact that the vortices are sitting in a periodic substrate. In the absence of a substrate the onset of diffusion is more gradual and a well defined melting temperature does not appear. Experimentally, the melting of the interstitial vortices occurs at or just below $`T_c`$  so that $`T_mT_c`$. In Fig. 10 we show a series of $`V`$ vs $`f_{DC}`$ curves for $`A=0.35`$ for $`T/T_m=`$ $`0.0`$, $`0.07`$, $`0.28`$, $`0.64`$, $`0.75`$, $`0.87`$, and $`0.96`$, from bottom to top. The fractional steps wash out very quickly with temperature at about $`T/T_m=0.1`$. The main step is visible all the way up to $`T_m`$ but is more cusp like for $`T/T_m>0.64`$. This washing out of the step with temperature is also consistent with the Shapiro step experiments near $`T_c`$ which found only a cusp feature . Experimentally it would still be possible to measure the increase of the step width with AC amplitude by taking the derivative of the $`IV`$ curves, which shows a dip at each side of the step. Since most transport experiments are performed close to $`T_c`$, it will be very difficult to observe the fractional steps. The widths of the steps can be increased by increasing $`A`$, allowing fractional steps to be visible for $`T>0.1T_m`$. In practice, however, there will be a limit to how large $`A`$ can be made before the vortices at the pinning sites begin to depin.

## IV Transverse Phase Locking For Triangular Pinning

Next we study the transverse phase locking for the case where the pinning is triangular rather than square. We use the same system size and pinning density as for the square pinning array, but every other pinning row is shifted in the $`x`$-direction by half a lattice constant. Here the motion of an interstitial vortex for zero AC amplitude

has a periodic transverse velocity component, unlike the square case where for zero AC amplitude the motion is strictly 1D. In Fig. 11 we show $`V`$ vs $`f_{DC}`$ for different values of $`A`$. We find phase locking similar to that seen in the square pinning case. The step widths all increase monotonically with AC drive, rather than oscillating as expected for Shapiro steps. There are some noticeable differences from the square case, in that the depinning force is only weakly affected by increases in $`A`$. In addition the main step is wider than that observed for square pinning at the same AC amplitudes. We also find that the step widths increase linearly with AC amplitude, rather the quadratically as in the square pinning array.

In Fig. 12(a) we show the vortex motion above depinning for $`A=0.0`$, indicating that the vortex moves with a periodic transverse component even in the absence of an AC drive. In Fig. 12(b-f) we show the vortex orbits along various steps for $`A=0.25`$. In Fig. 12(b), at $`f_{DC}=0.07`$, a stable periodic orbit forms where the vortex goes through an interesting loop that alternates from up to down. This loop feature appears for stable orbits where $`f_{DC}<0.1`$. For the step near $`f_{DC}=0.11`$ \[Fig. 12(c)\], corresponding to the second largest step, a stable orbit forms similar to that in Fig. 12(b) but the loop feature is absent. On the small step at $`f_{DC}=0.1375`$ \[Fig. 12(d)\], a stable orbit forms which is not symmetrical in the $`y`$-direction. A vortex trajectory for a non-step region at $`f_{DC}=0.2`$ \[Fig. 12(e)\] has similar characteristics as the trajectories for the non-step regions for square pinning. In Fig. 12(f) we show the vortex orbit along the main step, $`f_{DC}=0.3`$, where a well defined sinusoidal orbit can be seen as the transverse frequency generated by the DC motion matches one-to-one with the frequency of

the transverse AC drive. For higher DC drives, there are additional stable orbits along the steps as in the case for the square pinning. Similar noise signals for step and non-step regions as seen in Fig. 4 are also observed for the triangular pinning case.

In Fig. 13 we show the scaling of the width $`\mathrm{\Delta }`$ of the main step for triangular pinning with AC drive $`A`$, along with the width $`\mathrm{\Delta }`$ of the main step for square pinning for comparison. Here the magnitude of $`\mathrm{\Delta }`$ for the triangular pinning is greater than that for the square pinning for all AC values considered here. The width of the triangular steps also scales linearly with AC amplitude rather than quadratically as in the square pinning, so that for smaller AC amplitudes the difference in the magnitudes of $`\mathrm{\Delta }`$ becomes greater. This result suggests that the detection of transverse phase locking would be easiest for a system with triangular rather than square pinning arrays.

In Fig. 14 we show the scaling of the step width $`\mathrm{\Delta }`$ with AC amplitude $`A`$ for the fractional step at $`f_{DC}=0.11`$ (the orbit shown in Fig. 12(c)), along with the width of the main step and the change in value of the depinning force. Here the fractional step scales linearly with AC amplitude, just as the main step does. The difference in the depinning force at different values of $`A`$, which is much smaller than the widths of the main step and the fractional step, does not scale linearly with $`A`$ but increases quadratically or even faster.

We can compare our results for transverse phase locking in the triangular pinning to the results for random pinning. For random pinning  the vortices in the highly driven phase regain partial triangular ordering and a washboard signal appears. Unlike the square pinning case, in random pinning in the absence of a transverse AC drive the vortices do not move in strictly 1D channels but have a transverse velocity component centered around zero which can also have a washboard signal. For the random pinning it was found that the width of the transverse phase locking step oscillates as a function of AC amplitude, which is different from the results for the triangular or square pinning. For low AC amplitudes, however, the increase in the step width is linear with AC amplitude in the random pinning case. Using a simple model it was shown that any form of disorder that induces a transverse temporal order in the absence of an AC drive gives rise to a linear dependence of the step width for small AC amplitudes . This is similar to the case of the triangular pinning, where there is transverse temporal order present in the zero AC drive limit. In addition, for the case of random pinning, the vortex trajectories along the step show stable sinusoidal orbits in a similar manner to the stable orbits found in the periodic pinning cases.

## V Conclusion

We have investigated the transverse phase locking in square and triangular pinning arrays. For square arrays we find a series of fractional steps which correspond to stable vortex orbits. Along the non-step regions there are no stable periodic vortex orbits. All the fractional steps increase monotonically in width with AC amplitude. Along the steps a narrow band noise signal is present, while in the non-step regions the narrow band signal is broadened and large low frequency noise power appears. We have shown explicitly that the width of several of the fractional steps increases as the the square of the AC amplitude. The velocity vs drive near the depinning threshold and the main step scale with $`\beta =1/2`$.

We find that the phase locking is most pronounced for interstitial vortex filling fractions at which the interstitial vortices form a symmetrical ground state, such as at $`B/B_\varphi =1.25,1.5`$ and $`2.0`$. For filling fractions near these commensurate fillings, a partial phase locking occurs where certain regions of the sample have stable phase locked orbits while other regions are unstable and the vortex velocity in these regions increases linearly. For filling fractions where the ground states are disordered the phase locking is absent. With a finite temperature, the fractional steps appear only at low temperatures and are washed out at higher temperatures, while the presence of the main step can be detected up to the melting temperature $`T_m`$.

For triangular pinning arrays, where moving interstitial vortices have a periodic transverse motion even in the absence of a transverse AC drive, we again observe a transverse phase locking with all the step widths increasing monotonically with AC amplitude rather than oscillating as in the case of Shapiro steps. The depinning force is only weakly increased by increasing the AC amplitude. The main step for the triangular case is larger than that observed in the square pinning case, and its width increases linearly with AC amplitude, unlike the quadratic increase observed for the square case. The linear increase in the step widths with AC amplitude is similar to the results for transverse phase locking with random disorder with small AC drives, where there is also an ordered transverse motion in the absence of an AC drive. We show that the width of the fractional steps also increase linearly with AC amplitude, but the difference in the depinning force increases quadratically.

Our predictions should be testable for superconductors with periodic pinning arrays where only one flux line is captured per pinning site. The steps can best be observed in samples where Shapiro steps have already been observed with longitudinal DC and AC drives at filling fractions where a symmetrical vortex configuration occurs, such as $`B/B_\varphi =17/16`$, $`5/4`$, $`3/2`$, or $`2`$. For experiments performed near $`T_c`$, it is unlikely that the fractional steps can be observed; however, the main steps should be visible. A particular advantage of transverse phase locking steps over longitudinal (or Shapiro) steps is that the step width can be made arbitrarily large by increasing the AC amplitude or lowering the AC frequency. Our results should also be relevant for vortex motion in Josephson-junction arrays at commensurate fillings.

Acknowledgments: We thank D. Domínguez, N. Grønbech-Jensen, A. Kolton, V. Moshchalkov, D. Stroud, and L. Van Look for useful discussions. This work was supported by the US Department of Energy under contract W-7405-ENG-36.