text:
# High Spatial Resolution Mid-infrared Observations of Three Seyfert Galaxies

## 1 Introduction

There is much evidence (e.g., Krolik 1999) that nonthermal processes power the nuclei, i.e., the central sources, of Seyfert galaxies. In the mid-infrared, however, thermal radiation from heated dust contributes an additional emission component. This mid-infrared radiation, at characteristic temperatures of a few hundred kelvin, is believed to emerge from a disk-like structure that is heated by the inner accretion disk (Sanders et al. 1989). A question remains of the geometry of the central, presumably mainly nonthermal, source and how much of the surrounding dust is heated by the central source and how much is heated by starbursts which are located throughout the nucleus and thus form an extended cloud. It is difficult to resolve this dust component in direct observations using currently available telescopes because the angular extent of the region containing the emitting dust is small. Quantitatively, the Seyfert galaxies considered here have observed luminosities of 1010-1011L. If the dust is at a temperature of $``$300 K — the temperature at which blackbody emission peaks at 12.5 $`\mu `$m —and is heated by a central source, it is predicted to have a separation from the central source 2 - 20 pc or $`<`$ 0′′.1 depending on the composition of the dust. If, on the other hand, the dust were heated by hot stars in starburst regions, the extent of the region containing dust clouds could be much larger.

The diffraction limited full width at half maximum (FWHM) of the Keck 10-m Telescope — the largest aperture telescope presently available — is 0′′.31 at 12.5 $`\mu `$m. Deconvolution techniques and proper sampling of a nearby point source function (PSF) can improve this resolution significantly. For example, deconvolution applied to the image of NGC 1068 — the nearest and brightest Seyfert galaxy — have successfully resolved structure with a size of 0′′.1 ($``$7 pc at the distance of NGC 1068, Bock et al. 2000).

In this paper we present observations of the nuclei of three nearby, infrared bright Seyfert galaxies — NGC 1275, NGC 4151, and NGC 7469 — taken with the Keck 10-m Telescope. These observations, when deconvolved, approach the angular resolution necessary to distinguish thermal from nonthermal emission and, in one case, resolve the structure of the central source. In the mid-infrared, all three galaxies are known to contain a strong point-like nucleus and the observations were designed, and optimized, to probe this source. As a result, they are not optimized to delineate low surface brightness, extended emission. They are thus complementary to observations such as those by Radomski et al. (2003) of NGC 4151. We take the Hubble constant to be 75 km s-1Mpc-1.

## 2 Sample

The characteristics of the three Seyfert galaxies (Seyfert 1943), taken from the literature, are shown in Table 1. All are known to be bright in the mid-infrared (IRAS Point Source Catalog 1989; Miles et al. 1996). NGC 1275 and NGC 4151 are both compact in the mid-infrared. Miles et al. have shown, from observations at $``$10 $`\mu `$m with 0′′.7 resolution, that NGC 7469 consists of a compact nucleus surrounded by a lumpy ring $``$3′′ in diameter. The ring is presumably a region of active star formation. Miles et al. found that the nucleus plus ring structure in NGC 7469 accounts for $``$60% of the total 12 $`\mu `$m flux density of the galaxy measured by IRAS.

## 3 Observations

The observations were made using the Keck Long Wavelength Spectrometer (LWS; Jones & Puetter 1993) in the imaging mode at 12.5 $`\mu `$m ($`\mathrm{\Delta }\lambda `$= 1.16 $`\mu `$m) at the f/25 forward Cassegrain focus of the Keck I 10-m Telescope. The LWS Si:As array has 128 $`\times `$ 128 pixels. The scale is 0′′.08 pixel-1 giving a field of $``$10$`{}_{}{}^{\prime \prime }\times `$ 10′′. A log of the observations is given in Table 2.

In the mid-infrared there were no stars in the field from which to estimate point source profiles. Each measurement of the object size was therefore accompanied by multiple observations of one of the three bright stars, identified as the PSF in Table 2, which was nearby to the object and of comparable infrared brightness. The observations were made in runs of 27 – 54 s by interleaving imaging of the galaxy with observations of the PSF star; the number of comparisons between the galaxy and the PSF star is also indicated in Table 2. The telescope was guided during the observations of the objects using a visual guider on stars near the galaxies.

For these observations, objects and PSF calibrators were imaged in a chop-nod scheme where the throw of the secondary chopper was equal to the amplitude of the (telescope) nod. The amplitude of the secondary chopper was set to 5′′ – 15′′ at a frequency $``$5 Hz.

All three nights were photometric and in all cases the airmass was between 1.0 and 1.4. Each sequence to measure the object size was preceded by an observation of one of the three bright stars — $`\alpha `$ Tau=HR 1457, $`\alpha `$ Boo=HR 5340, or $`\beta `$ Peg=HR 8775 — in order to obtain a photometric calibration.

## 4 Data Reduction and Analysis

### 4.1 Photometry

Photometric measurements were made from the raw images before any deconvolution algorithms (described below) were applied. The photometry of the galaxies was based on the following magnitudes (in the Vega based system) — HR 1457: \[12.5 $`\mu `$m\]=-3.07 mag, HR 5340: $`[12.5\mu `$m\]=-3.15 mag, and HR 8775: $`[12.5\mu `$m\]=-2.55 mag. The “sky” was taken as an annulus 4$`\stackrel{}{\mathrm{.}}`$6 to 6$`\stackrel{}{\mathrm{.}}`$0 in diameter. The conversion from magnitude to flux density followed the prescription given in the Explanatory Supplement to the IRAS Catalog (Beichman et al. 1985); i.e., at 12.5 $`\mu `$m, 0.0 mag was taken as 26.15 Jy.

### 4.2 General Image Data Reduction

The images were processed following a similar procedure to that described in Bock et al. (2000), briefly summarized as follows. Images were obtained by shifting and combining sub-images after every chop pair, using only the central image, to eliminate any possibility of non-reproducibility in the telescope nod, i.e., potential differences caused by guiding the telescope were obviated in software by finding the centroid of each object image in each nod. Both end chop images did not always fall onto the array, so were not used in the data reduction. The results when the end images were included, however, were found to be appreciably the same.

For the PSFs, which were not guided, each chop-nod combination consisted of four separate images, which were combined to produce three images: an image near the middle of the array and two images near the edge or off the array. The middle image was thus seen twice, first in the positive beam of the first chop pair, and then in the negative beam of the second chop pair. The objects were guided, and eight sub-images were combined into the three images; the middle image was thus seen four times.The object and PSF calibrator images were sub-pixel shifted, rotated, and co-added using a routine that relies on sub-Nyquist sampling of the images. The objects were rotated from array coordinates to equatorial north coordinates during the combination. The PSF calibrators were rotated through a sequence of rotations in array coordinates matched to the object. The rotations were designed to best match the orientation of the characteristic six-pointed diffraction pattern as it appeared on the array for both object and PSF calibrator. This pattern is fixed with respect to the telescope and thus in paralactic angle and therefore rotates in the field as the telescope tracks. The resulting PSF was thus not rotated to equatorial north, but was rotated through a superposition of rotation angles corresponding to the orientation of the diffraction pattern of the object sequence as it appeared on the array. The rotations for each PSF/object combination were always matched pair-wise, producing a matched PSF for every subset of object images under consideration.

### 4.3 FWHM Analysis

Because the combined images of the objects were all very nearly point-like and closely resembled the combined PSF calibrator images, and variations in the PSF, not statistical noise, dominated the overall uncertainty in deconvolving these images, it was necessary to evaluate the stability of the PSF by comparing the raw object and PSF calibrator images before proceeding to image deconvolution. An estimate of the intrinsic source size of the object was obtained by assuming the sizes add in quadrature,

$`\theta _{INT}^{}{}_{}{}^{2}=\theta _{OBJ}^{}{}_{}{}^{2}\theta _{PSF}^{}{}_{}{}^{2}`$,

where $`\theta _{INT}`$ is the intrinsic source size of the object, $`\theta _{PSF}`$ is the observed PSF size, and $`\theta _{OBJ}`$ is the observed object size all expressed as FWHM. An estimate of the uncertainty can be gotten from $`\delta \theta _{INT}^{}{}_{}{}^{2}`$, the dispersion of $`\theta _{INT}^{}{}_{}{}^{2}`$ in the population of pairs. We will take the quantity $`\theta _{UL}=[<\theta _{INT}^{}{}_{}{}^{2}>+\delta \theta _{INT}^{}{}_{}{}^{2}]^{1/2}`$ as the estimate of the upper limit on the source size. It should be noted that this definition takes account of the fact that the objects and PSFs are often comparable in size, and the fact that the size of the object image on occasion exceeds that of the PSF.

### 4.4 Deconvolution

Each combined object image was deconvolved by a matched PSF calibrator image using the Richardson-Lucy procedure (Richardson 1972; Lucy 1974) as encoded in the IDL maximum likelihood routine (Varosi & Landsman 1993). The images were interpolated to a four times finer pixel scale before deconvolution. The reconvolved image was rebinned and subtracted from the original image and compared with the estimated statistical uncertainty per pixel. It was found, however, that the uncertainties were largest where the image was brightest, and were always greater than the statistical uncertainties even after iterating the deconvolution routine a large number of times. The non-uniform uncertainties from deconvolution must be compared to a combination of statistical and PSF uncertainties in the region of image of interest. Therefore we chose to stop iterating the deconvolution routine when the residuals were consistent with the measurement uncertainties in the raw images.

## 5 Results

NGC 7469 showed a statistically significant detection of structure while the images of NGC 1275 and NGC 4151, combining all the data, both showed a single point-like compact nucleus . The FWHM of various angular quantities, obtained before deconvolution was applied, are given in Table 3. Although the uncertainties on the $`<\theta _{INT}^2>^{1/2}`$ are small, we claim, to be conservative, only to have measured a limiting size for the nuclei of NGC 1275 and NGC 4151 of $`\theta _{UL}`$. It should be noted that the PSF during the NGC 7469 observations was noticeably more compact and significantly more stable than for any of the other objects. Deconvolution applied to the combined NGC 7469 image confirmed the presence of clear structure.

In the case of NGC 1275 and NGC 7469, the data were extensive enough to provide several sets of the object and the PSF calibrator. The images were paired to both match the object and PSF calibrator images in airmass and to minimize the time interval between the two images. Unfortunately, the NGC 4151 data were quite limited, resulting in only two image pairs. In order to obtain a meaningful limit, the data were simply sub-divided into five image pairs. This is less than ideal since the pairs all have a common time interval, and the matching in airmass was poor. Further sub-dividing the NGC 4151 images was possible, but gave a nearly identical result. From the present data we conclude that 0′′.16 is an upper limit to the size of the nucleus. In comparison, Neugebauer et al. (1990), using slit scans and a single element detector, determined the size of the nucleus of NGC 4151 at 11.2 $`\mu `$m to be 0′′.16 $`\pm `$ 0′′.04. The deconvolved images of NGC 4151 and NGC 1275 do not show structure and these sources appear to be unresolved in these observations. As stated above, the deconvolution algorithm did not improve the size limit beyond that set by the variations in the PSF and we choose to quote the upper limit of the sizes as the size estimated in the analysis of the FWHM in the raw images.

As expected from the work of Miles et al. (1996), the raw images of NGC 7469 show a bright compact nucleus surrounded by a fainter ring. The flux from the ring, which has a diameter of $``$3′′ and is well separated from the nucleus, is comparable to that of the nucleus. The surface brightness from the ring is less than 5% of the brightness of the central source in the raw image. \[Note that because of the small chop amplitude of these observations we cannot make a reliable flux measurement in an adequately large beam to determine this ratio.\]

Photometric results are included in the summary of Table 4. With 20 % uncertainties, the flux densities of the compact sources of NGC 1275 and NGC 4151 reported here at 12.5 $`\mu `$m agree with those measured at 12 $`\mu `$m by IRAS in a 4.5 $`\times `$ 0.8 beam. These results are consistent with the results on NGC 4151 reported by Radomski et al. (2003), who find that 73 % of the 10 $`\mu `$m broad band flux is in an unresolved component. We will assume that the luminosity and color as measured by IRAS applies to the compact sources. The flux density of the compact nucleus of NGC 7469 inside a 1$`\stackrel{}{\mathrm{.}}`$4 diameter beam is 40 % of the flux density at 12.0 $`\mu `$m measured in the much larger IRAS beam.

### 5.1 Deconvolution of NGC 7469

The deconvolved image of NGC 7469, after combining all the data, shows significant structure both in the ring and in the nucleus. The deconvolution routine reaches statistical noise on the 3′′ diameter ring in about ten iterations, due to its low surface brightness. Maximum likelihood does not perform well on images like this requiring large dynamic range, and we note that although many more iterations are needed to resolve the nucleus, they are of no further use in resolving the ring. The deconvolution of the image with ten iterations, thus emphasizing the 3′′ diameter ring, is shown in Figure 1. For purposes of deconvolving the nucleus, the ring does not affect the images of the central source. This assertion was tested by deconvolving a 3$`\stackrel{}{\mathrm{.}}`$$`\times `$ 3$`\stackrel{}{\mathrm{.}}`$8 image of NGC 7469 containing the ring, and by deconvolving a 1$`\stackrel{}{\mathrm{.}}`$$`\times `$ 1$`\stackrel{}{\mathrm{.}}`$6 image of NGC 7469 excluding the ring. Both deconvolutions produce essentially identical images of the inner 1$`{}_{}{}^{\prime \prime }\times 1^{\prime \prime }`$ nucleus. In the following we will often designate the central nucleus of NGC 7469 simply as “NGC 7469”.

The deconvolved image of the NGC 7469 nucleus with 1000 iterations of the deconvolution algorithm is presented in Figure 2 and shows extended structure that is statistically significant based on the estimate of PSF variations. The robustness of the deconvolution technique was tested by sub-dividing the observations of the PSF into subgroups and treating successive subgroups as objects. This test is described in Appendix A. The reproducibility of the extension shown in Figure 2 was likewise tested by sub-dividing the data into independent object and PSF calibrator images and deconvolving the corresponding image pairs, as shown in the Appendix B. The statistical significance can be evaluated quantitatively by comparing the quality of fit of the hypothesized elongated source to that obtained under the assumption of a simple point source. This procedure is described in Appendix C with the result that the elongation is estimated to be significant at the 4.3-$`\sigma `$ confidence level.

To summarize the appendices, we conclude the extension of ($`<`$0$`\stackrel{}{\mathrm{.}}`$4) $`\times `$ 0$`\stackrel{}{\mathrm{.}}`$08 seen in NGC 7469 is real for the following reasons:

1. Structure in NGC 7469 was detected statistically based on the FWHM of the object and the PSF.

2. Deconvolved PSF images, synthesized from the observations for matching in time and airmass, are point-like and more compact than the deconvolved images of NGC 7469.

3. Deconvolved images using independent pairs of NGC 7469 and the PSF produce consistent results with a northwest-southeast extension of similar size and direction.

### 5.2 Extended Emission

Although the observations were not optimized to study extended emission, the quantity and quality of the data on NGC 1275 and NGC 7469 were sufficiently great to attempt an examination of the extended emission using the technique exploited successfully by Radomski et al. (2003) in the study of NGC 4151. These high resolution images, combined with a precise knowledge of the PSF, allow us to best probe extended emission close to the nucleus.

In Figure 3, two sets of differences – the galaxy minus the PSF and the galaxy minus a fraction of the PSF –are given. There is evidence in NGC 7469 for emission in the northwest–southeast direction — i.e., the orientation of the structure in the deconvolved image — extending almost 1′′ while in NGC 1275 there is an indication of extended emission in the north-south direction. The putative extended emission is only about three times brighter than the PSF residuals, so we advise caution in interpreting these images, especially at the radius of the first Airy maximum (radius = 0.4′′). Although faint, the extended emission reproduces itself if the images are constructed using independent subsets of the data from the first half and second half of the observations, in spite of the fact that the change in paralactic angle rotates the apparent diffraction pattern by 10 $`\mathrm{°}`$. In NGC 7469, the structure noted here is also seen in Figure 1.

## 6 Discussion

The observations at 12.5 $`\mu `$m with the 10-m Keck Telescope yield a resolved structure in NGC 7469, but give only limits for the size of NGC 1275 and NGC 4151; see Table 4. At the distances of the galaxies, the limits on the angular sizes correspond to linear size limits of 10 pc in NGC 4151 and 27 pc in NGC 1275. Although an asymmetry is obvious in NGC 1068 at a comparable physical scale, none is obvious in NGC 4151 in our images. The physical temperatures implied for dust grains heated by a central source of luminosity are between 300 K and 450 K depending on the sizes of the individual grains and their composition. This is substantially hotter than the observed color temperature, based on IRAS 12 $`\mu `$m and 25 $`\mu `$m data, of 192 K and 201 K for NGC 1275 and NGC 4151 respectively. As shown in Table 4, the limits on the sizes in these sources are not small enough to require optically thick emission, but the lower bound on $`\tau _{12.5\mu m}`$ is 0.25 in NGC 1275 and 0.06 in NGC 4151 if the emission is thermal emission at the color temperature determined from the IRAS 12 $`\mu `$m to 25 $`\mu `$m flux density ratio. This calculation, of course, requires the reasonable assumption that the source structure is essentially the same at 25 and 12 $`\mu `$m.

In the case of NGC 7469, the resolved structure provides strong direct evidence for the dust thermal emission model in the nucleus as well as a limit of $`\tau _{12.5\mu m}>`$ 0.4 within this structure (Table 4). The elongated structure is suggestive of a disk geometry. The “standard model”of an active galactic nucleus (AGN) (e.g., Krolik 1999; Sanders et al. 1989) suggests that this disk is heated by the central AGN. The fact, noted above, that graphite and silicate dust grains at distances of $``$15 pc from the central source, i.e. at the linear radius corresponding to the angular size of the structure, are heated to temperatures of 300 - 450 K, well above the observed color temperature, is consistent with this picture. The derived brightness temperature, again well below 300 K, is only a limit since the source is unresolved in one dimension.

An interesting question is what is the relationship of the infrared nucleus to the nuclear structure of the galaxy observed in the radio. Observations by Thean et al. (2001) at 1.6-GHz (18-cm) at $``$ 0$`\stackrel{}{\mathrm{.}}`$3 resolution show structure whose orientation is very similar to that seen in the mid-infrared structure. A similar feature is also seen at 4.9-GHz (6-cm) in the VLA image of Wilson et al. (1991). Although again at a much lower spatial resolution (2′′) than that seen here in the mid-infrared, high resolution maps in CO (Meixner et al. 1990) show emission, interpreted as a rotating molecular disk, whose orientation also is similar to that seen in the mid-infrared structure.

VLBI observations of Lonsdale, Lonsdale & Smith (2003) at 18-cm show that the nucleus of NGC 7469 is a nearly symmetric linear triple source aligned in the east-west direction with the outer components separated by 0$`\stackrel{}{\mathrm{.}}`$15, and the central source nearly centered between these two components. The agreement between the centroids of the radio peak, as described by Lonsdale et al. , and the mid-infrared peaks is excellent, the two peaks coinciding to within 0$`\stackrel{}{\mathrm{.}}`$040 if the 3′′ ring structures are assumed to be spatially coincident. The uncertainty in the determination of the astrometric coincidence of the radio and infrared nuclei is determined by how well the structures in the 3′′ ring can be aligned.

We have compared the mid-infrared image of Figure 2 with the VLBI image of Lonsdale et al.(2003) in Figure 4. The position angle of the elongated structure at 12.5 $`\mu `$m is 135°, $``$45° off of the position angle of the linear radio structure. To register the images in Figure 4 we assumed that the middle radio peak coincides with the AGN nucleus, and is centered on the mid-infrared peak. Then the two outer radio components are symmetrically placed with respect to the mid-infrared disk, and is suggestive of radio jets emerging out of the plane of the dust disk. If the galaxy is assumed to be circular, its apparent size on the sky – 1.1’ $`\times `$ 1.5’ (de Vaucouleurs et al. 1991) – implies it is at an inclination angle of $``$45°.

The surface luminosity densities of the three compact nuclei, as listed in Table 4, are on the order of a few times 1014 L Kpc-2. Soifer et al. (2000) found that, of seven Ultraluminous Infrared Galaxies (ULIRGs) observed in the mid-infrared with the Keck 10-m Telescope, only Markarian 231, a known Seyfert 1 galaxy, has a luminosity density as high as 1014 L Kpc-2. The difference between starbursts and AGNs is heightened by an examination of the surface brightnesses of infrared luminous galactic nuclei listed by Soifer et al. , Soifer et al. (2001), and Evans et al. (2003). Galactic HII regions have surface brightnesses which range from 2 $`\times `$ 1011 to 2 $`\times `$ 1012 LKpc-2, those of the starburst galaxies studied range from 2 $`\times `$ 1011 to 1013 LKpc-2, while those of the ULIRGS, with the exception of Markarian 231, range from 2 $`\times `$ 1012 to 6 $`\times `$ 1013 LKpc-2. Gorjian, Turner & Beck (2001) have reported the infrared surface brightness of a “super” star cluster of $``$ 3 $`\times `$ 1014 L Kpc-2 over 1 – 2 pc in the nucleus of NGC 5253. While the peak surface brightness in the “super” star cluster matches that of the Seyfert galaxies the physical extent of the “super” star cluster is smaller, consistent with the 1 – 2 orders of magnitude lower luminosities compared to the Seyfert nuclei. Thus it is seen that the infrared surface brightness of those galaxies with evidence of AGN activity is significantly higher than that of those dominated by star formation and we can conclude that star formation provides only a negligible fraction on the luminosity of the nuclei of the Seyfert galaxies measured here. These results argue that the infrared surface brightness, in concert with the infrared luminosity, can be used as a diagnostic measure for the presence of starbursts or AGNs. These observations can be used to measure the relative contributions of AGN and starbursts to the energy budgets of infrared luminous galaxies provided that the AGN is not heavily obscured at the observed wavelength.

## 7 Summary

High spatial resolution observations at 12.5 $`\mu `$m were made with the 10-meter diameter Keck Telescope in order to study the central nuclei of three Seyfert galaxies.

1) The observations, when deconvolved, resolve a linear structure in the Seyfert galaxy NGC 7469 that is ($`<`$13) $`\times `$ 26 pc. The orientation of the infrared emission, when compared to CO and radio observations, supports the conjecture that the mid - infrared emission traces a disk.

2) The observations do not detect structure in the nuclei of NGC 1275 (FWHM $`<`$ 28 pc) or of NGC 4151 (FWHM $`<`$ 10 pc). The former limits are within a factor of two of being able to resolve a significant contribution to the total emission by thermal emission of centrally heated dust grains.

3) There are indications of underlying extended emission of about 300 pc in NGC 7469 and NGC 1275.

4) The luminosity surface brightness derived for the three Seyfert galaxies are a few times 1014 LKpc-2, significantly higher than that associated with starbursts.

We thank the staff of the Keck Observatory for their assistance in making these observations possible. The W. M. Keck Observatory is operated as a scientific partnership between the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. It was made possible by the generous financial support of the W. M. Keck Foundation. We extend special thanks to those of Hawaiian ancestry on whose sacred mountain we are privileged to be guests. Without their generous hospitality, none of the observations presented herein would have been possible. We thank Carol Lonsdale for supplying the VLBI.image of NGC 7469. T.S. and E.E. were supported by grants from the NSF and NASA. B.T.S. and L.A. are supported by the SIRTF Science Center at Caltech. SIRTF is carried out at J.P.L., operated by Caltech under an agreement with NASA.  APPENDICES
A - PSF STABILITY  The stability of the PSF accompanying NGC 7469 was estimated by deconvolving a synthesized image of the PSF calibrator by a complementary set of independent images of the PSF calibrator. It is important to accurately match the two PSF images in time and airmass. When deconvolving the NGC 7469 image, the average airmass of the object differs only by 0.005 from the average airmass of the PSF, and the average time interval between object and PSF sub-images is eight minutes. Subgroups to serve as the object and PSF were synthesized using the fact that each of the four sub-images of the PSF actually consists of four independent chop-nod sets; in Table A1 these are called a, b, c, and d. The composition of several synthesized PSFs is shown in Table A1. Each pair of PSFs are independent, and were chosen with variable matching of the PSF pair in time and airmass. Deconvolved images resulting from the use of one subgroup as PSF and one as object are illustrated in Figure 5. Close matching in time and airmass (e.g., PSF\_C//PSF\_D) produce point-like deconvolved images. Poor matching in airmass (e.g., PSF\_12//PSF\_34) results in significant extensions. Matching in airmass seems to be somewhat more important than matching in time, as the PSF\_13//PSF\_24 and PSF\_12//PSF\_34 pairs produce the most extended deconvolved images. For example, the pair PSF\_14//PSF\_23 has the same time interval as the pair PSF\_12//PSF\_34 but much better matching in airmass. The most representative matching for purposes of deconvolving the object are PSF\_C//PSF\_D.  B - NGC 7469  The image of NGC 7469 was deconvolved by combining all of the object and PSF data. The object and PSF images are very closely matched in airmass and time interval (see Table A2). The data set was divided into independent images of the object and PSF to assess the stability of the northwest-southeast extension in the deconvolved image of NGC 7469. Splitting the data into two independent images of NGC 7469 and the PSF, N\_12//PSF\_12 and N\_34//PSF\_34, results in two images which both show a similar northwest-southeast extension; see Figure 6. This agreement is remarkable given that the image pairs are completely independent. As shown in Figure 6, further sub-dividing the data into four sub-images and PSFs results in four individual images which all show the northwest-southeast extension. Furthermore, all images also show two smaller isolated structures located about 0′′.5 to the northwest and southwest. A smaller structure located $``$0′′.3 to the east is seen in the two image pairs, but not in each of the four independent sub-images.  C - Statistical Significance of the Structure  The statistical significance of the elongated structure apparent in the deconvolved images of NGC 7469 was evaluated by comparing the quality of fit of the hypothesized elongated source to that obtained under the assumption of a simple point source. This comparison is based on the residuals with respect to the observed data for the two cases, taking into account our knowledge of the properties of the measurement noise. The residual images are shown in panels (e) and (f) of Figure 7, along with the sequence of images used to generate them, consisting of: (a) observed image, (b) PSF, (c) “model” image representing the central source obtained by deconvolution, and (d) the model image convolved with the PSF. The residual image (e) was obtained by subtracting (d) from (a), and is to be compared with (f), which represents the residual image corresponding to the best point-source fit, i.e, the residual image obtained by subtracting (b) from (a). The comparison of (e) and (f) shows that the residuals are noticeably lower for the extended source than for the point source. The relative probabilities of the two source-geometry hypotheses can be calculated by comparing the conditional probabilities of the extended source and a point source given the data (represented by a vector y whose components represent the pixel values of the observed image) by $`P(\mathrm{ext}|𝐲)`$ and $`P(\mathrm{point}|𝐲)`$. A simple application of Bayes’ rule then yields:
$$\mathrm{ln}P(\mathrm{ext}|𝐲)/P(\mathrm{point}|𝐲)=\frac{1}{2}\underset{i}{}[(y_ic_i^{\mathrm{point}})^2(y_ic_i^{\mathrm{ext}})^2]/\sigma _i^2$$
(1)
where $`c_i^{\mathrm{point}}`$ and $`c^{\mathrm{ext}}`$ represent the point-source and extended-source model images convolved with the PSF, respectively, and $`\sigma _i`$ represents the measurement noise in the $`i`$th pixel. The summation is over all statistically independent pixels in the image. In deriving this equation a flat prior was assumed, representing equal a priori probabilities of extended and point sources.

The measurement noise is assumed to be a Gaussian random process with variance:

$$\sigma _i^2=\sigma _{\mathrm{bg}}^2+f^2V_i$$
(2)
where $`\sigma _{bg}`$ is the standard deviation of the background (obtained using pixels far from the source itself), $`f`$ is the source strength in data units, and $`V_i`$ is the “variance map” of the PSF, representing the normalized value of PSF uncertainty as a function of position. The PSF variance map was obtained using the deconvolution residuals of four independent sub-images of the PSF calibrator and smoothing the result using a $`5\times 5`$-pixel boxcar averaging window. The PSF uncertainty term in the above expression represents a good approximation in the case of a marginally resolved source.

Based on this procedure, a probability ratio $`P(\mathrm{point}|𝐲)/P(\mathrm{ext}|𝐲)=1.5\times 10^5`$ was obtained. After normalizing the sum of the two probabilities to unity, this result implies that the probability that the data were produced by a point source was of this same order, which corresponds to a $`4.3\sigma `$ (or greater) deviation of a Gaussian random variable. We can therefore equate the significance of our result to a $`4.3\sigma `$ detection of source elongation.


format: text

text:
# Atomic-scale perspective on the origin of attractive step interactions on Si(113)

## Abstract

Recent experiments have shown that steps on Si(113) surfaces self-organize into bunches due to a competition between long-range repulsive and short-range attractive interactions. Using empirical and tight-binding interatomic potentials, we investigate the physical origin of the short-range attraction, and report the formation and interaction energies of steps. We find that the short-range attraction between steps is due to the annihilation of force monopoles at their edges as they combine to form bunches. Our results for the strengths of the attractive interactions are consistent with the values determined from experimental studies on kinetics of faceting.

Self-assembly of steps is an important way to achieve surface patterning at length scales where the usual lithographic techniques are not applicable. On certain vicinal semiconductor surfaces, the equilibrium bunching of steps leads to the formation of periodic hill-and-valley patterns (or grooves), that can serve as templates for growing regular arrays of one-dimensional (1D) nanoscale devices like quantum wires. Experimental work in recent years mochriePRL ; mochrieGrooveSS ; mochrieSS ; sudohPRL ; zippingLetter  has revealed that such groove structures are formed on vicinal Si(113), where the flat (113) surface coexists with step-bunched regions over a wide range of temperatures and misorientation angles.

In a series of experiments, Mochrie and coworkers mochriePRL ; mochrieGrooveSS ; mochrieSS  have studied the orientational phase diagrams for stepped Si(113) surfaces miscut towards several crystallographic orientations. In all the cases, their studies strongly suggest that, at low temperatures, the steps gather in bunches as a result of attractive step-step interactions. Independent experimental work on the fluctuations of the single and multiple-height steps and the kinetics of faceting transitions by Sudoh and coworkers sudohPRL ; zippingLetter  have confirmed the presence of short-range attractive interactions between steps. Furthermore, their analysis of the step-zipping process has provided quantitative estimates for the magnitude of the attractive interactions zippingLetter . The presence of such interactions has been assumed in theoretical models that show the orientational phase diagram mochriePRL  can be explained in terms of competing short-range attractive and long-range repulsive step interactions lassig ; vbsPRL . While much progress has been made in understanding the self assembly of steps on Si(113) surfaces, a microscopic analysis of the step interactions has not been attempted and the physical origin of the short-range attraction has not been convincingly explained.

In this letter, we investigate the interactions between steps on Si(113), with the goal of understanding the origin of the effective short-range attraction between steps. Describing each step as a pair of equal and opposite force monopoles located at the step edges, we find that the short-range attraction between steps is due to the cancellation of force monopoles that belong to adjacent steps in a step bunch. This observation is confirmed by investigating the scaling of the step formation energies and interaction strengths with the step height. Noting that the structure of the steps is determined by the direction of the miscut, we consider here Si(113) surfaces miscut towards $`[\overline{3}\overline{3}2]`$, so that our estimates for the magnitudes of step interactions can be directly compared with the recent experiments of Sudoh et al. zippingLetter .

The slab geometry that shows the Si(113) terraces and types of steps studied in the present work is given in Fig. 1. Steps run along the $`[1\overline{1}0]`$ direction and consist of (114) nanofacets zippingLetter , with a ($`2\times 1`$) reconstruction erwin114 . For the (113) surface we consider both the (3$`\times `$1) and the (3$`\times `$2) structures dabr  that are stabilized by interstitial atoms. We distinguish two possible structures for the (3$`\times `$2) case, depending on the position of the interstitial rows on the terraces with respect to the step. When the interstitials on the lower terrace are away from the step or close to the step, we call them the (3$`\times `$2)a and the (3$`\times `$2)c structures, respectively (refer to Fig. 1).

The dimension of the simulation cell in the $`x`$-direction, $`L_x`$, is determined by the number $`k`$ of Si(113) unit cells and the number $`n`$ of Si(114) unit cells that constitute the step: $`L_x=n\mathrm{\Delta }+ka_x`$, where $`a_x=12.73`$Å and $`\mathrm{\Delta }=16.21`$Å are, respectively, the periodic lengths of the (113) and the (114) surfaces along $`[\overline{3}\overline{3}2]`$ (Fig. 1). The periodic length in the $`[1\overline{1}0]`$ direction $`L_y`$ is the smallest length compatible with the periodicity of both the Si(113) and the Si(114) structures, $`L_y=23.04`$Å. The thickness of the slab $`L_z`$ was chosen to be $``$150 Å. To simulate the vicinal geometry, we use shifted periodic boundary conditions farid  in which the shift along  is determined by the step height $`nh`$, where $`h=1.64`$Å is the height of a single step.

The energies of the slabs in Fig. 1 are computed using empirical and tight-binding (TB) models for the atomic interactions abinitio . Empirical calculations are carried out with the Tersoff potential T2 tersoff2  as it reproduces the structure, surface stress and relative energy of the Si(113) and Si(114) reconstructions reasonably well. While the tight-binding methods are more accurate, they are usually impractical for ledge energy calculations since they require slabs with several thousands of atoms to capture the long-range interactions between steps. To circumvent this problem, we have developed a multi-scale approach where the tight-binding model is coupled to an empirical potential as described in the following paragraph.

The key idea of the coupling scheme is to handle the long-range elastic fields produced by steps using an empirical potential whose elastic constants are very close to their TB counterparts; this would ensure that the elastic fields in the hybrid scheme would not be very different from the fields obtained in the case where all the atoms were treated using TB. We use a charge self-consistent TB method, whose superior transferability has been studied in detail in czwJPCM . This method is coupled with the Tersoff (T3) potential tersoff3  using the scheme shown in Fig. 2. We divide the simulation cell into three regions: (1) a thin surface region where the complex bonding topologies in the vicinity of the step-edge are treated with the TB method, (2) a thick ”bulk-like” region where all the atoms with tetragonal bonding environments are treated with the T3 method and (3) a padding zone which is used to communicate atomic displacements between the first two regions. The atoms in the simulation cell are relaxed using the following procedure (refer to Fig. 2): (1) relax the thin slab (Zone 1, Fig. 2(a)) using TB, while keeping the padding zone fixed; (2) update the atomic coordinates of Zone 1 in Fig. 2(b) with the values obtained after the TB relaxation of Zone 1 in Fig. 2(a); (3) relax the thick slab and the padding zone (Zone 2 and P.Z. in Fig. 2(b)) using the T3 potential while keeping the atoms in Zone 1 fixed; (4) update the coordinates of the P.Z. atoms in Fig. 2(a) with the values obtained after the T3 relaxation of P.Z. in Fig. 2(b); (5) repeat steps (1) through (4) until energy convergence is achieved elsewhere .

Once the total energy $`E`$ of the slab is calculated, step interactions can be obtained from the ledge energy farid

$$\lambda (L_x)(ENe_b\gamma _{113}L_xL_y)/L_y,$$
(1)
where $`e_b`$ is the bulk cohesion energy per atom ebulk , $`N`$ is the number of atoms in the cell in Fig. 1, and $`\gamma _{113}`$ is the surface energy of Si(113).

Using the fact that the steps interact through long-range dipolar fields Mdipoles , we have analyzed the short-range interactions by fitting the ledge energies in Eq. (1) to the form vbsPRL

$$\lambda (L_x=n\mathrm{\Delta }+ka_x)=\mathrm{\Lambda }_n+G_n\frac{\pi ^2\mathrm{\Delta }^2}{6L_x^2}U_n\delta _{k0},$$
(2)
where $`\mathrm{\Lambda }_n`$ is the formation energy per unit length of a step of height $`nh`$, $`G_n`$ is the strength of the long-range dipolar interactions, and $`\delta _{k0}`$ denotes the Kronecker delta symbol, where $`k`$ is the number of (113) unit cells on the terrace (Fig. 1). The energy offset $`U_n`$ that appears in Eq. (2) describes the effective interactions at the smallest possible spacing between the steps, $`n\mathrm{\Delta }`$ (note that the (113) terraces are absent in this case). A positive (negative) value of $`U_n`$ indicates the presence of a short-range attraction (repulsion) between steps. Our results for different terrace and step structures are discussed below.

We first focus attention on the vicinal surfaces made of single-height steps and illustrate the general features of the step interactions by considering the energetics of the (3$`\times `$2)c steps obtained from the T2 potential. In Fig. 3, we plot the ledge energy as function $`1/L_x^2`$ ($`L_x`$ is the step separation) and provide the fitting parameters $`\mathrm{\Lambda }_1`$, $`G_1`$ and $`U_1`$ defined in Eq. (2). If the dipolar description of step interactions is assumed to hold over the entire range of step separations (as implied by Eq. (2)), Fig. 3 shows the presence of a short-range attraction ($`U_1`$), since the point corresponding to $`L_x=\mathrm{\Delta }`$ lies below the linear fit for the ledge energy. The physical reason for the origin of this attraction will be discussed next.

We start by noting that the Si(113) and Si(114) facets are atomically compatible, in the sense that they can meet to form edges parallel to \[1$`\overline{1}`$0\] without any changes in their individual bonding topologies. Subsequently, the elastic fields of a (114) step (or nanofacet) can be expressed in terms of a dipole which consists of a pair of equal and opposite force monopoles ($`F_+,F_{}`$)Mgrooves  at the edges of the step, as shown in inset (a) of Fig. 3. While such steps are expected to show inverse-square repulsive interactions at large separations, we find that this behavior holds for separations as small as 29Å (this point is marked in Fig. 3). When the steps are closer than 29Å, they lose their individuality forming a smooth Si(114) surface: at a separation of $`L_x=\mathrm{\Delta }`$ there is a cancellation of the monopoles at the step-edges as illustrated in inset (b) of Fig. 3. The absence of dipoles with mutual repulsive interactions is the reason for the presence of an effective short-range attraction between steps. Further evidence for this argument is obtained by computing the interactions of multiple-height steps. We will focus on such calculations after we discuss the interaction parameters for other types of steps given in Fig. 1.

Using both the T2 potential and the TB-T3 method, we find the existence of short-range attractive interactions between steps for all the reconstructed surfaces considered here. However, the magnitudes of these interactions depend strongly on the structures of the terraces and steps (refer to Table 1). Since the TB method allows for the tilting of the dimers $`d_1`$ and $`d_2`$ on the step (marked in Fig. 1), we have reported the results for the cases where the tilting patterns are in-phase and out-of-phase restrict . The magnitude of the attractive interaction $`U_1`$ is $`60`$meV/Å and $`10`$meV/Å for the TB-T3 method and T2 potential, respectively. This discrepancy between the values of $`U_1`$ computed with TB-T3 and with T2 is due to the differences in step formation energies obtained with the two methods, as well as to the sensitivity of $`U_1`$ with respect to changes in $`(\gamma _{114}^0\gamma _{113})`$ (refer to Fig. 3). For both T2 and TB-T3, the repulsion strength $`G_1`$ is of the order of 10meV/Å, while the formation energy $`\mathrm{\Lambda }_1`$ is of the order of 100meV/Å (see Table 1). In order to compare our results with the experimental ones zippingLetter , we note that Sudoh et al. zippingLetter  report the quantity $`\epsilon \mathrm{\Lambda }_1\mathrm{\Delta }(\gamma _{114}^0\gamma _{113})`$, i.e. the difference between the step formation energy $`\mathrm{\Lambda }_1`$ and the excess energy of the Si(114) facet per length of the step. Table 1 shows that the values for $`\epsilon `$ obtained with the TB-T3 method are of the same order of magnitude as the estimates of Sudoh et al. zippingLetter . Closer agreement could be achieved if the long-range interactions are included in the analysis of kinetics of step-zipping reported in zippingLetter . Also, future experiments aimed at determining the structure of the step ($`a`$ or $`c`$) could verify the trends given in Table 1.

We will now focus on the formation energy and interactions of steps with height $`nh`$ ($`n`$-bunches) for $`n=2,3`$. Because of the larger size of the simulation cells, we have used only the T2 potential tersoff2  for this analysis. In the zero temperature limit, the formation energy of an $`n`$-bunch can be expressed as the sum of the formation and interaction energies of the steps within the bunch vbsPRL :

$$\mathrm{\Lambda }_n=n\mathrm{\Lambda }_1U_1(n1)+G_1\underset{i=1}{\overset{n1}{}}\frac{ni}{i^2}.$$
(3)
Due to the presence of attractive interactions between single-height steps, the formation energy of an $`n`$-height step is expected to show sub-linear dependence of $`n`$, $`\mathrm{\Lambda }_n<n\mathrm{\Lambda }_1`$ vbsPRL . This prediction is indeed borne out by the atomistic results given in Table 2. Using the parameters $`\mathrm{\Lambda }_1,G_1,U_1`$ from Table 1 for the ($`3\times `$2)c structure, we find that the bunch energies predicted by Eq. (3) are within 1–2 meV/Å of the values for $`\mathrm{\Lambda }_n`$ obtained in the present simulations.

Turning our attention to the step-step interactions, we find that the repulsive interactions $`G_n`$ show a nearly quadratic dependence of the step height $`nh`$, while the strength of the attractive interaction $`U_n`$ shows a weak dependence weakUn  on $`n`$ (refer to Table 2). The quadratic behavior of $`G_n`$ provides strong evidence to the idea regarding the cancellation of the force monopoles associated with the adjacent steps in an $`n`$-bunch. When such cancellation takes place, the only remaining monopoles are associated with the ends of the bunch. Consequently, each $`n`$-bunch behaves as a force-dipole with a moment that is proportional to $`nh`$. Identical $`n`$-bunches interact as dipoles Mdipoles , with a repulsion strength determined by the square of the dipolar moment of one $`n`$-bunch, which results in the $`n^2`$ dependence of $`G_n`$.

In conclusion, we have shown that the effective short-range attraction experienced by steps on Si(113) can be explained by the cancellation of force monopoles associated with adjacent steps in a step bunch. While numerical estimates for the interaction parameters depend on the interatomic potential used, our physical reasoning for the origin of attractive interaction is robust and is consistent with experiments. A key point in understanding how the atomic structure can lead to step-step attractive interactions is the atomic compatibility of the two stable surface orientations, Si(113) and Si(114): when nanofacets with these orientations intersect forming a step edge, no bond breaking or rebonding occurs. Recent experimental work takeguchi  indicates that no rebonding occurs at the step-edges formed between the Si(113) and Si(337) facets; our calculations show the presence of attractive interactions between steps in this case, further lending support to our physical picture. The difference between the type of steps considered here, and the most common steps on Si(001) surfaces, such as the SB and the DB steps chadi  is that when steps on Si(113) merge they form a new facet (e.g. Si(114)), while such atomically smooth facets are not observed in the latter case due to significant rebonding at the step-edges chadi . Therefore, a simple cancellation of force monopoles is not expected for the SB and the DB steps on Si(001), and hence they do not show attractive interactions farid . Further theoretical and experimental studies aimed at elucidating the dependence of the step-step interactions of the azimuthal orientation will shed more light on the self-assembly of steps on Si(113) miscut in an arbitrary direction.

Support from the MRSEC at Brown University (DMR-0079964), NSF (CMS-0093714, CMS-0210095), and the Salomon Research Award from the Graduate School at Brown University is gratefully acknowledged. Ames Laboratory is operated by the U.S. Dept. of Energy and by the Iowa State University (contract W-7405-Eng-82).